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BACKGROUND OF THE INVENTION 

As the capabilities of microprocessors have increased, hardware modules 
dedicated to IEEE floating-point division have became common in microprocessors. 
Appendix A.-Table 1 lists the latencies (i.e., number of cycles required to complete an 
instruction) and restart times (i.e., number of cycles that elapse until a functional 
module can engage in a new independent computation) of floating-point division 
modules in various microprocessors. It is easy to see that floating point division is 
rather slow compared to addition and multiplication. The relevant IEEE standard is 
IEEE Standard for Binary Floating-Point Arithmetic, ANSI/IEEE 754- 1 985 

When division is performed by numerical iteration using Newton's well- 
known method, it is necessary to perform two multiplication operations, one of which 
is dependent on the result of the other, thereby requiring that the multiplication 
operations be performed sequentially. This requirement limits the speed and 
efficiency of division operations using Newton's method. In contrast, the well-known 
algorithm of Goldschmidt, which also requires two multiplication operations, relies on 
two multiplication operations which are independent of one another, and which can 
therefore be performed simultaneously to improve the speed and efficiency of the 
division. In terms of processor architecture, Goldschmidt's division algorithm is more 
amenable to pipelined and parallel implementations. 

Unfortunately, however, there is currently in the prior art no satisfactory 
measure of the error when employing the Goldschmidt algorithm, and without a good 
measure of error, the inaccuracies of the intermediate iterative computations 
accumulate and cause the computed result to drift away from the correct quotient. 
That is, implementations of Goldschmidt's algorithm are not self-correcting. The lack 
of a general and simple error analysis of Goldschmidt's division algorithm has 
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deterred many designers from implementing GoIdschmidt ! s algorithm. Thus, most 
implementations of multiplicative division methods have been based on Newton's 
method in spite of the longer latency due to dependent multiplications in each 
iteration. 

Those who have utilized Goldschmidt's algorithm have had to keep careful 
track of accumulated and propagated error terms during intermediate computations. 
Recent implementations of Goldschmidt's division algorithm still rely on an error 
analysis that over-estimates the accumulated error. Such over-estimates lead to correct 
results but require a costly full-precision multiplier circuit that wastes hardware 
resources and causes unnecessary delay, because the multiplier and the initial lookup 
table are too large. 

The lack of an accurate error estimator therefore discourages the use of 
Goldschmidt's division algorithm, and prevents an efficient realization of the pipeline 
advantages of Goldschmidt's algorithm when implemented. This results in increased 
hardware complexity, power consumption, processing time, and latency for division 
operations. 

There is thus a widely recognized need for, and it would be highly 
advantageous to have, a compact, accurate, and efficient error estimator for a 
Goldschmidt division algorithm conforming to IEEE formatting and rounding 
standards. This goal is met by the present invention. 
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SUMMARY OF THE INVENTION 

The present invention is of an apparatus and method conforming to IEEE 
formatting and rounding standards, which efficiently and accurately estimates error 
when using Goldschmidt's algorithm. This allows efficient use of pipelining to 
increase the speed of floating-point division operations. 

In an embodiment of the present invention, multiplication is performed by a 
Booth recoded multiplier that can be fed by: 

(a) a redundant Booth recoded operand and a non-redundant binary operand; 

(b) two redundant carry-save operands; or 

(c) two non-redundant binary operands. 

In another embodiment of the present invention, IEEE rounding is performed 
by a novel "dewpoint" rounding technique that implements dewpoint rounding with 
back multiplication. Performing back multiplication with an estimated dewpoint 
allows the use of a half-size multiplier instead of a full-size multiplier in estimating 
Goldschmidt algorithm error. 

Accordingly, yet another embodiment of the present invention utilizes a half- 
size multiplier in performing Goldschmidt's algorithm, yielding IEEE-correct 
rounding while offering the advantages of: 

(a) reduced hardware; 

(b) improved pipeline organization; 

(c) fewer clock cycles; 

(d) shorter clock cycles; 

(e) reduced latency; 

(f) increased throughput; and 

(g) lower power consumption. 
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Therefore, according to the present invention there is provided a method for 
IEEE-rounding a computed quotient in a processor, the computed quotient 
corresponding to an exact quotient which equals a dividend divided by a divisor, the 
method including: (a) determining an error range of the exact quotient; (b) 
determining a first candidate number and a second candidate number from the error 
range; (c) associating the first candidate number with a first rounding interval 
containing numbers that are IEEE-rounded to the first candidate number; (d) 
associating the second candidate number with a second rounding interval containing 
numbers that are IEEE-rounded to the second candidate number; (e) computing the 
dewpoint number, which separates the first rounding interval from the second 
rounding interval; (f) back-multiplying the dewpoint number by multiplying the 
dewpoint number by the divisor; and (g) comparing the back-multiplied dewpoint 
number against the dividend to determine whether the first candidate number 
represents the IEEE-rounded computed quotient, or whether the second candidate 
number represents the IEEE-rounded computed quotient. 

Furthermore, according to the present invention there is provided a Booth 
multiplier for computing the product of a first operand and a second operand, 
including: (a) a first stage operative to preparing the first operand and the second 
operand for the addition of partial products, and operative to recoding the second 
operand in Booth radix-8 digits, and operative to generating partial products; (b) a 
second stage having an adder tree operative to compressing the partial products; and 
(c) a third stage having an adder operative to compressing the carry-save 
representation of the product to a binary representation. 
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DESCRIPTION OF THE PREFERRED EMBODIMENTS 

The principles and operation of an apparatus and method according to the 
present invention may be understood with reference to the accompanying description 
and the material in the appendices, which disclose the detailed mathematical 
principles, circuit diagrams, examples, and other information to completely explain 
implementing the invention. 
"Dewpoint" Rounding 

The present invention discloses a novel rounding procedure for IEEE floating 
point division (see Appendix D-Section 4), which is herein referred to as "dewpoint" 
rounding — so-called because of the analogy to the meteorological temperature below 
which condensation of moisture occurs, and above which condensation of moisture 
does not occur. The procedure relies on an error range of the quotient that allows for 
only two candidate numbers for the final IEEE rounded result. Each candidate number 
is associated with a rounding interval, which is simply the set of numbers that are 
IEEE-rounded to the candidate number. The dewpoint is defined to be the number 
separating the two rounding intervals. The rounding decision is obtained by 
comparing the dewpoint against the exact quotient by applying back-multiplication. 
This comparison determines which of the candidate numbers is the correct IEEE- 
rounded result. Appendix D-Section 4 discloses the details of a unified dewpoint 
rounding procedure for all IEEE rounding modes, thereby eliminating the need for 
rounding tables. 

The novel dewpoint rounding of the present invention represents an 
improvement over current prior-art approaches for implementing IEEE rounding in 
division, which first compute a "rounding representative" of the exact quotient (i.e., a 
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number that belongs to the same rounding interval to which the exact quotient 
belongs), and then round the rounding representative. 

An optimized implementation of dewpoint rounding and back multiplication is 
disclosed in detail in Appendix B-Section 7.6 and in Appendix D-Section 5.3. 

To compute the dewpoint, use a rounding injection and a dewpoint 
displacement constant (as detailed in Appendix D-Section 4.2). The rounding 
injection is added to the computed quotient, which is then truncated. The dewpoint 
displacement constant is then added to the truncated computed quotient to obtain the 
dewpoint. 

All the intermediate results mentioned here (the computed quotient, the 
truncated computed quotient, and the dewpoint) are also represented in redundant 
representation. This means that each of the additions mentioned above can be 
computed in constant time (i.e., they do not require a carry-propagate adder with a 
logarithmic delay). This enables a reduction of the four IEEE rounding modes to a 
single rounding mode, so that each separate mode need not be dealt with separately. 

Furthermore, the test to determine which of the candidate numbers represents 
the proper rounding of the computed quotient involves evaluating whether the 
quantity (b*dewpoint - a) is zero, positive, or negative. It is noted that the dewpoint is 
very close to the exact quotient a/b 9 so that the absolute value of this quantity is very 
small, and the sign (or zero) of this quantity can be determined by the least-significant 
few bits. Thus, the hardware used to determine the sign (or zero) can be fed by a small 
subset of the least-significant bits. 

Moreover, in yet another embodiment of the present invention, the back- 
multiplication is split into two half-size multiplication operations that can be 
performed in consecutive clock cycles using the same multiplier (refer to cycles 8 and 
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9 in Appendix B-Table 3). The first part of the back-multiplication is done with an 
estimated dewpoint (as shown in Appendix B-Section 7.6). The estimated dewpoint is 
computed from the computed quotient of the previous iteration. An apparatus 
according to this embodiment of the present invention utilizes a half-size multiplier 
for the dewpoint back-multiplication, thereby reducing hardware requirements. 

Multiplier Hardware Optimization 

Addition trees in multipliers are not amenable to pipelining. Short clock cycles 
are therefore not achievable at reasonable cost if the addition tree has too many rows. 
Booth radix-8 recoding reduces the number of rows in an addition tree from n to 
(«+l)/3. 

Booth radix-8 multipliers are usually implemented using a 3-stage pipeline, as 
follows: 

precompute the 3x multiple of the first operand of the multiplier and 
recode the second operand; 

an addition tree that computes a carry-save representation of the product; 
and 

final carry-propagate addition. 

Goldschmidt's algorithm performs only two multiplications per iteration. 
Hence running Goldschmidt's algorithm on a 3-stage pipeline creates unutilized 
cycles ("bubbles") in the pipeline. These bubbles increase the latency and reduce the 
throughput. Certain prior-art processors attempt to utilize these bubbles (and increase 
throughput) by allowing other multiplication operations to be executed during such 
bubbles. 

The present invention discloses a Booth radix-8 multiplier that allows for both 
operands to be either in nonredundant representation or carry-save representation. 

01494376X11-01 



' 9 

Booth multipliers with one operand in redundant carry-save representation are known 
in the prior art, but Booth multipliers with both operands in redundant representation 
conceptually is a novel feature of the present invention, which reduces the 3-stage 
pipeline to a 2-stage pipeline for all but the last iteration of the algorithm. 

The Booth-8 multiplier design that supports operands in redundant 
representation is not symmetric, in the sense that the first operand and second operand 
of the multiplier are processed differently during the first pipeline stage. During the 
first pipeline stage, operands represented as carry-save numbers are processed as 
follows: 

(a) The first operand is compressed and the 3x multiple thereof is computed. 
This requires two adders. 

To reduce hardware requirements, an embodiment of the present invention 
employs the adder from the third pipeline stage for compressing the first 
operand. In this embodiment, the compression of the first operand appears 
in Appendix D-Table 2 as an operation that takes place in the third 
pipeline stage. 

(b) The second operand can be partially compressed from carry-save 
representation before being fed to the Booth recoder. In Appendix D- 
"Implementation of the dewpoint computation" a recoding method is 
detailed. 

A method for determining the Booth recoding of the dewpoint correction term 
is detailed in Appendix B-Section 7. This method is based on a bound on the value of 
the dewpoint correction term, which determines the most significant digit position / 
involved in the computation (i = 24 in Appendix B-Figure 3). 
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A first Booth recoded operand of the dewpoint correction term is computed 
modulo 2'', which has either the value of the dewpoint correction term or the value of 
the dewpoint correction term plus 2"'. A second Booth recoded operand is computed in 
the same manner, minus 2"'. Only the most significant Booth recoded digit of the 
second Booth recoded operand needs to be computed. The other digits are the same as 
in the first Booth recoded operand. 

In parallel with the above computations, a signal is computed that indicates 
whether the first Booth recoded operand represents the dewpoint correction term plus 
2"'. If the signal is a zero, the first Booth recoded operand represents the dewpoint 
correction constant, and is chosen as the Booth recoded operand. If the signal is a one, 
the Booth rcoded operand represents the dewpoint correction constant plus 2"' and the 
second Booth recoded operand is chosen as the Booth recoded operand. For example, 
2"' = 2' 24 in Appendix B-Figure 3. 

In an embodiment of the present invention, a Booth recoded multiplier can be 
fed by either non-redundant binary operands or by redundant carry-save operands. 
When applied to a Booth radix-8 multiplier, this enables reducing the feedback 
latency to two cycles. The prior art features only Booth multipliers with one operand 
in redundant carry-save representation or signed-digit representation. 

The organization of a Booth multiplier according to this embodiment of the 
present invention has the following stages, as detailed in Appendix D-Section 5.1. 

Stage 1. The two operands of the multiplier are prepared for the addition of 
the partial products in the second stage. The second operand is 
recoded in Booth radix-8 digits and the partial products are 
generated. If the second operand is given in carry-save 
representation, a partial compressor prepares the second operand for 
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the input of a conventional Booth recoder. The recoding can accept 
either a binary string or a carry-save encoded digit string. 
The first operand is processed as follows: The 3x multiple of the 
operand is computed using an adder. The first operand can be 
represented in either binary or redundant carry-save representation. 
For the case where the first operand is encoded as a carry-save digit 
string, the computation of the 3x multiple is preceded by a 4:2 adder 
that computes a carry-save encoding of the 3x multiple. This carry- 
save encoded digit string is compressed to a binary number by the 
adder. For the case where the first operand is encoded as a carry- 
save digit string, the binary representation of the operand is also 
computed by a binary adder. The binary adder from the third 
pipeline stage can be used for this purpose if available for carry-save 
feedback operands. This can save an adder in the first pipeline stage. 

Stage 2. In the second stage, the partial products are compressed by an adder 
tree. In addition to the partial products, an additional row can be 
dedicated for an additive input. 

Stage 3. The third stage contains an adder to compress the carry-save 
representation of the product to a binary representation. This adder 
can be shared with the first pipeline stage. 

Appendix D-Section 6.1 details how a full size multiplier is used in a floating 
point divider. Appendix B-Section 6 details how a half-sized multiplier is used. 

While the invention has been described with respect to a limited number of 
embodiments, it will be appreciated that many variations, modifications and other 
applications of the invention may be made. 
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Appendix A 

Apparatus for Pipelined Division with IEEE 
Rounding 

(This appendix was part of the provisional patent application 60/421,998) 

Abstract 

We propose optimized pipelined implementations for Goldschmidt's division algorithm with IEEE 
rounding based on Booth Radix-8 multiplication. The considered optimizations for the quotient 
approximation are based on a careful general analysis of tight error bounds for the implementa- 
tion and are accompanied by the utilization of redundant representations, partial compressions, 
injection-based rounding and rectangular multipliers for the internal computations. For the imple- 
mentation of IEEE rounding, we introduce the concept of dewpoint rounding, that allows efficient 
implementation and reduced requirements for the quotient approximation. 

On this basis we propose the implementation of different versions of Goldschmidt's division 
algorithm with different pipeline depths. None of these implementations requires a full-sized mul- 
tiplier at any stage of the computations. In this way we reduce latency, cost, and enable increased 
throughput at a reasonable cost. We suggest a full range of pipelining depth: On one extreme is 
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a 3-stage pipeline with a restart time that simply equals the latency minus the number of pipeline 
stages. On the other extreme is a fully pipelined design. 

A.l Introduction and Summary 

As the number of transistors in microprocessors increased, a hardware module dedicated to IEEE 
floating point division became common in microprocessors. Table A.l lists the latencies (i.e., 
number of cycles required to complete an instruction) and restart times (i.e., number of cycles that 
elapse till a functional module can engage in a new independent computation) of floating point 
division modules in various microprocessors. One can readily see that floating point division is 
rather slow compared to addition and multiplication. In addition, designs based on multiplicative 
division contain a full precision multiplier, a rather costly circuit. 

This paper presents implementations of Goldschmidt's floating point division algorithm [16] 
with significantly reduced latency and cost. These implementations suggest a reasonable cost- 
pipelining tradeoff that allows for a full range of restart times. On one extreme is a three-stage 
pipeline with a restart time that simply equals the latency minus the number of pipeline stages. On 
the other extreme is a fully pipelined design (i.e. restart time of one cycle). 

Both extremes are quite practical. The largest multiplier used in any pipeline stage is half- 
sized, and hence the implementations are amenable to short clock periods. The cost of the designs 
is also reasonable. For example, in IEEE double precision, the fully pipelined design contains five 
rectangular half-sized multipliers (i.e., size 62 x 30) and two quarter-sized multipliers (i.e., size 
62 x 15. The number of partial products associated with these multipliers is roughly 3.6 times the 
number of partial products associated with a full-sized multiplier (i.e. 56 x 56). On the other hand, 
the three-stage pipeline design contains only a single half sized multiplier. 

We refer the reader to Appendix C and the references for a short review of division algorithms. 
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Related work. Goldschmidt [16] developed a floating point division algorithm in which each 
iteration requires two independent multiplications. The algorithm maintains three numbers TV* 
(which converges to the quotient), (which converges two 1), and F { = 2 — £)* (called the 
scaling factor). The algorithm converges quadratically. Loosely speaking, this means that the 
number of bits needed to represent Fi doubles in each iteration. This leads to the motivation of a 
fully pipelined design using a sequence of multipliers (one or two per iteration), where the second 
operand's length doubles in this sequence. 

Probably the main hurdle in implementing Goldschmidt's algorithm is analyzing the errors 
due to imprecise intermediate computations. In Appendix C a general error analysis of a version 
of Goldschmidt's algorithm is presented. The advantage of the error analysis is demonstrated by 
analyzing a micro-architecture similar to that in [28] and showing that a smaller multiplier could be 
used. An error analysis similar to that in Appendix C enables one to realize the idea of a sequence 
of multipliers (one per iteration), where only the last multiplier is "full-sized". We further carry 
an idea of Goldschmidt utilizing the fact that the scaling factor tends to 1 so that even the last 
multiplier is "half-sized". 

Techniques. The first step in this paper is an extension of the error analysis in Appendix C. One 
reason for having to extend the analysis is the usage of redundant representations (e.g. carry-save) 
for intermediate results. Redundant representation incurs two problems: (i) an increase of relative 
errors due to truncation (which is easily handled by the analysis in Appendix C), and (ii) inability 
to determine if an intermediate result is less than (or greater than) 1. However, our main motivation 
for extending the error analysis is to guarantee that the scaling factor F { (which tends to 1) never 
falls below 1. For this purpose we introduce an assumption which we call saturated rounding of F{. 
Loosely speaking, the fact that the scaling factor Fi tends to one implies that the most significant 
half of the binary representation of Fi is either 1.00 ... or 0.11 ... [16]. This observation implies 
that a half-sized multiplier suffices. 

The second step in this paper is a new rounding procedure for IEEE floating point division (see 
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Sec. A.4). We refer to this rounding procedure as dewpoint rounding. The procedure relies on an 
estimate of the quotient that allows for only two rounded results. A rounding interval corresponds 
to each candidate rounded result. The point that separates these rounding intervals corresponding 
to the candidate rounded results is called the dewpoint. The rounding decision is obtained by 
comparing the dewpoint and the exact quotient. This method differs previous approaches that (i) 
compute a rounding representative of the exact quotient (i.e. a number that belongs to the same 
rounding interval that the exact quotient belongs to), and (ii) round the representative. We present 
a unified dewpoint rounding procedure for all rounding modes that avoids rounding tables. 

We employ additional techniques, among them: (i) A Booth recoded multiplier that can be fed 
by either non-redundant binary operands or by redundant carry-save operands [7]. This technique 
when applied to a Booth radix-8 multiplier enables reducing the feedback latency to two cycles. 

(ii) Injection based rounding is used to implement directed rounding of intermediate results [9]. 

(iii) Back-multiplication (i.e comparison between the dewpoint and the exact quotient) requires 
full precision multiplication. Since we only have a half-sized multiplier, two passes through the 
multiplier are required. We rearrange the sequence of multiplications so that one of these passes 
takes place during what is otherwise a bubble in the pipeline. Hence, the latency associated with 
back multiplication is only one cycle. 

Organization. In Section A.2 notation is introduced, division is defined, and Goldschmidt's al- 
gorithm is presented. An error analysis of Goldschmidt's algorithm appears in Section A.3. This 
analysis builds on the analysis that appears in Appendix C. In Section A.4 dewpoint rounding is 
presented. Section A.5 lists various techniques employed to save cost and delay in a hardware im- 
plementation. Section A.6 presents a hardware design, the pipeline organization and scheduling, 
and evaluates design options. 
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Table A.l: Latencies(restart-times) of floating-point operations in current commercial micropro- 
cessors compared to the proposed division implementations. 

A.2 Preliminaries 

Notation Let XiX i+1 • • * xj £ {0, 1}* denote a binary string. We often denote this string also by 
x[i : j]. We also sometimes refer to Xi as x[i]. Since we deal with fractions (mostly in the binade 
[1,2)), the weight associated with the bit Xi is 2~\ Namely, a fraction is represented by the binary 
digit string X0.X1X2 • • • x p -\. 

Let a = <Ji<Ti_|_i • • • dj 6 {0, 1, 2}* denote a carry-save encoded digit string (i.e. <7i 6 {0, 1,2}). 
A carry-save encoded digit string x[i : j] is represented by two binary vectors xc[i : j] and xs[i : j]. 
Each carry-save digit x[£] satisfies x[£] = xc[£] +xs[£]. We refer to xc (resp. xs) as the carry-string 
(resp. sum-string) of X. 

A carry-save encoded digit string a[z : j] may represent two values: an unsigned value and a 
signed value (often called a two's complement value). The unsigned value represented by a[i : j] is 
denoted by \a[i : j]\ and equals \a[i : j}\ = Y^t^i a ^~ t - The two's complement value represented 
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by the carry-save representation a[i : j) is denoted by \a[i : j]\ two and equals \cr[i : j]\ two = -<n • 

Division Operation We consider the following task which captures the main difficulty in IEEE 
floating-point division. The dividend and the divisor are represented by p-bit binary strings A[0 : 
p - 1] and B[0 : p - 1], where \A\, \B\ G [1, 2). Our goal is to compute Q[0 : p - 1] that is the 
binary encoding of the rounded value of the normalized quotient. In other words, (i) Let A = A 
if \A\ > \B\, and A f = leftshift(A) if \A\ < \B\. This way, the quotient |A'|/|£| e [1, 2). (ii) let 
q = |j4'|/|i?|, (iii) round q (according to the IEEE standard). The binary string Q[0 : p — 1] should 
satisfy \Q\ = q. 

The precision determines the value of p as follows. In double precision operation p = 53, and 
in single precision p = 24. Four rounding modes are defined in the IEEE standard: round towards 
zero round towards zero (RZ), round to nearest even (RNE), round towards infinity (RI) and round 
towards minus infinity(RMI). 

Goldschmidt's division algorithm In Appendix C a variation of Goldschmidt's division algo- 
rithm based on directed rounding is presented and analyzed. The algorithm is listed below as 
Algorithm 1. The multiplications involved and their dependencies are depicted in figure A.l. 

Directed roundings are used for all intermediate calculations. For example, N- is obtained by 
rounding down the product N^ x • F[_ x . We denote by the relative error incurred when N-_ x • F{_ x 
is rounded down. Rounding down translates to the assumption that rc* > 0. Similarly, rounding 
down is used for computing F[ and rounding up is used for computing D\. Our error analysis 
requires that: (a) the operands are in the range A \ B satisfy A G [B z 2B) and B G [1, 2), so that 
A jB G [1, 2) (b) all the relative errors incurred by directed rounding be at most 1/4, and (c) we 
require that |e 0 | + 3do/2 + /o < 1/2. 
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Figure A.l: Schedule of the Iterations of Goldschmidt's division algorithm using approximate 
arithmetic and an initial approximation for 1/5. The numbers in circles indicate the sequence 
of the multiplications involved. For our implementation a bound on the relative error p(N-) of 
iteration i appears in each iteration. 



A.3 Extension of error analysis 

In this section we present a variation of the error analysis presented in Appendix C. The error 
analysis presented in Appendix C used a simplifying assumption called strict directed rounding 
(SD-rounding). Satisfying SD-rounding is easy if intermediate results are represented in non- 
redundant binary representation. In our division algorithm redundant representation is used for 
intermediate results (i.e. carry-save and borrow-save encoded digit strings). Hence we cannot rely 
on the SD-rounding assumption. 

The analysis uses the same requirements used before, namely: (a) the operands A 1 and B satisfy 
A* 6 [B, 2B) and B e [1, 2), (b) all the relative errors incurred by directed rounding are at most 
1/4, and (c) |e 0 | + 3d 0 /2 + / 0 < 1/2. 



24 



Algorithm 1 Goldschmidt-Approx-Divide(A', B) - Goldschmidt's division algorithm using ap- 
proximate arithmetic 

Require: 0 < n* + d t + /> < 1/4, for all i, and |e 0 | + 3d 0 /2 + /o < 1/2. 
1: Initialize: TV^ <- A 7 , D'^ <- B, Fi x <- 
2: for z = 0 to k do 



6: end for 

7: Return(Ar/) 

A.3.1 Rounding assumptions 

The strict-directed rounding assumption is defined as follows. 

Assumption 1 fSD-roundingJ Directed rounding up (resp. down) r is strict if x < 1 implies 
r(x) < 1 (Vesp. x > 1 implies r(x) > 1). 

Instead of assuming that SD-rounding holds, we use the following assumptions for iterations i > 0: 
Assumption 2 (saturated rounding) Fori > 0, 



Note that the definition of the relative error fi in the computation of F( is now ambiguous (i.e. 
suppose F[ = 1 and (1 — fi) • (2 — D^) < 1). To avoid ambiguity, define the relative error f[ as 
follows: 



Note that if D\ < 1, then > // > 0. This means that saturated rounding reduces the relative 
error in the computation of F[ when < 1. 

Assumption 3 (precise multiplication by 1) Every multiplication by 1 is precise. 



3: ^(l-O-iV^.i^. 
5: ^'^(l-/ i ).( 2 -^). 



J^-maxf!, (!-/«)• (2 -££)}. 
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The assumption that multiplication by 1 is always precise means that if F- = 1, then = £> t ' 
and jV/ +1 = N<. We say that the algorithm is stuck starting with iteration i if Dj = D[ and 
Nj = for every j > i. Note however, that F{ = 1 does not imply that the algorithm is stuck 
starting with iteration z. The algorithm may not be stuck if < 1 and eventually the relative error 
fj is small enough so that F( > 1. However, if Z>- > 1, then the algorithm is stuck starting with 
iteration i. 

In the following intermediate values obtained with saturated rounding and precise multiplica- 
tion by 1 are simply denoted by TV/, D\, and F{. 

A.3.2 Analysis 

With the following theorem we show that even with our current rounding assumptions, N[ approx- 
imating A 1 /B from below and the relative error p(N<) = A '#^** can be bounded from above by 
the same expression as in the error analysis in Appendix C. 

Theorem 4 For every i > 0, the relative error p(Nl) = A jrjjf* satisfies 



0< p{N[) < TTi+Si 



(A.1) 



where TTi is defined by 7Ti = 1 — (1 



ni)H 




l f: > 0 and where Si is defined by 




|e 0 | -h 3do/2 fori = 0 
+ fi-i otherwise 



for i > 0. 



Lemma 5 The following three statements hold: 

a) IfFU = (1 - /,_!) • (2 - DU) andDU e [1 - 5^,1 + *_J f then D\ > 1 - 

b) IfFU = (1 - / w ) • (2 - A'-i) *«i £>i < 1 + *■ 
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c) IfFU = 1 > (1 - fi-i) • (2 - DU) and Di. x e [1 - 6^, 1], fAen D< € [1 - 1]. 
Proof: Let tj = — 1. We first prove Parts a) and b): 

D\ = (l + *)-(l + ti-i)-(l-*_i)-(l-./i_i) 



= (1 + di) ■ (1 - t?_i) • (1 - /i-i) > (1 - d - fi-i) = 1 - <$* this implies a) 

>1 <1 <1 

< + this implies b) 



c) Since F(_ x = 1, precise multiplication by 1 implies that D\ = D'i_ l9 so all we need to show 
is that D'i > 1 - S € . Since 1 > (1 - f^ x ) • (2 - £> t '_i), we can write: 



A / >(l-/i-i)-(2-A , -i)'^-i 

^(i-z^o-a-ei) 



and Part c) follows. □ 

Proof of Theorem 4: In Appendix C Equation (A.l) is proven for iteration i whenever D[_ x e 
[1 - 6i- U 1 + Si-i] and F{_ x = (1 — • (2 — D'^). In fact, these two conditions are used in 
Appendix C to prove that 1 - Si < • F{_ x < 1, which implies Equation (A.l) for iteration i. 
We rely on the proof presented in Appendix C and deal here with the effect of saturated rounding. 
Let io and ip be defined as follows: 



i D = min{z > 0 : D\ > 1} 

i F = min{2 > 0 : (1 - /<) • (2 - < 1} 



Note that if D\ > 1, then (1 - /*) • (2 - D$ < 1, and hence, i F < i D . 
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We need to consider the 3 cases: (i.) i < i F \ (ii.) i F < i < <z>; and (iii.) i > i D . Note that case 
(ii.) only occurs if i F 

Case (L): One can directly verify that 1 - S 0 < D' Q < 1 + 6q. Since saturated rounding does not 
take place in iteration i = 0, it follows that Fq = (1 — / 0 ) • (2 - D' Q ). Hence the error analysis in 
Appendix C implies that Equation A.l holds for i = 1. 

For 0 < i - 1 < i F we have D'^ < 1 and F{_ x = (1 - /,_!) ■ (2 - D^) > 1 since 
i - 1 < min{i D ,i F ), By applying induction and Part a) of Lemma 5, it follows that 6 [1 — 
1]. Hence the proof in Appendix C applies, and Equation A.l holds for i. It follows that 
Equation A.l holds for every i € [1, i F ). 

Case (ii.): By applying induction and Parts a) and c) of Lemma 5, it follows that D[_ x 6 [1 — 
1], for every i - 1 < i D . If F{_ x = (1 - • (2 - D-^), then Equation (A.l) holds for L 
If F[_ x = 1 > (1 - /;_!) • (2 - D<_i), then by Part c) of Lemma 5, D\ <E [1 - 1]- By precise 
multiplication by 1 it follows that D[ = and hence 1 — Si < • < 1. Hence the error 

analysis in Appendix C is applicable, and Equation (A.l) holds for i < i D . 

Case (iii.) (i > io): Note that the algorithm is stuck starting from iteration %d- Since 7r J+ i > 7Tj, 
it suffices to prove that 0 < p{N' iD ) < 7r( ii?+1 ). 

We now prove that D iD e [1, 1 + d io ]. The lower bound follows from the definition of i D . 
Part b) of Lemma 5 proves that if = (1 - fi D -i) • (2 - D/^-i), then D io < 1 + di D . So 

all we need to prove is that F{ D _ X = (1 - fi D -i) • (2 - £>•,,_!). Since D iD > 1 it follows that 
if D iD -i < 1, then Ff D _ x ^ 1, hence = (1 - f iD _ x ) - (2 - Z^.j). The minimality of Z£> 

implies that if Ad-i > 1, then = 1. But saturated rounding is not used in iteration i = 0, so 
the upper bound on D[ D holds. 
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N' 

We now expand -^f as follows: 



Thus, 



_ A rr l-ll( 



and 



Since D' iD > 1, it follows that 



1 — 7i£ 



<i- ( i-„ ((B+1) ,n^ 



lJ? 1 



which gives the desired upper bound on p{N[ D ). To prove the lower bound we use D' iD < 1 + d iL 
with equation (A.2) as follows: 

1 — Tl£ 



^ e )>i-(i + 4»)IIi^ 



£=0 
(io-1) 



= *io > 0. 
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This completes the proof for the last of the three cases. □ 
Corollary 6 For every i > 0, F( € [1, 1 + 

Proof: If F( = 1, then the corollary obviously holds. We deal now with the case that F- = (1 — 
fi) • (2 - Di). In the proof of Theorem 4 we showed that D\ > 1 - 5*, for i > 0. Therefore, 
F/ = (1 — fi) • (2 — DJ) < 1 + <5i. The lower bound follows from Assumption 3 (saturated 
rounding). □ 



A.4 IEEE rounding of the quotient 
A.4.1 Background on IEEE rounding 

The IEEE Floating-Point Standard 754 [17] requires the implementation of all four IEEE round- 
ing modes for all basic arithmetic operations and for all precisions supported. IEEE rounding is 
supposed to be computed on the exact result of the operation (which is A'/B in our case). For 
most operations the exact result is not available or too expensive to be computed. It is therefore 
common practice to compute the IEEE rounding result by computing a rounding representative 
first and compute IEEE rounding on the representative that leads provably to the same result [10]. 

In particular for multiplicative division it is even difficult to find a representative quotient from 
an approximation of the quotient, and it is common practice to involve a back-multiplication by B 
of the approximation and a comparison with A 1 for this purpose [19, 34]. 

A different approach is suggested in [18], which avoids the back-multiplication, but involves 
computing roughly twice the accuracy of the target precision and making conclusions about round- 
ing from this representation. This approach is based on a proof that shows that, in the binary 
representation of the exact quotient, the length of runs of zeros (ones) must be limited. We will 
not follow this approach since it involves one additional iteration for multiplicative division algo- 
rithms, adds significant delay to the implementation, and requires very wide intermediate operand 
representations for the last iteration. 
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In the sign-magnitude representation of floating-point numbers, the four IEEE rounding modes 
can be reduced to three rounding modes RZ, RI, and RNU (round to nearest up) based on the 
numbers' sign [32]. 

Observation 7 The exact quotient A'/B cannot be a midpoint between two representable num- 
bers. Therefore RNU and RNE are equivalent rounding modes with respect to division. 

Based on the above reduction and observation we only need to consider the three rounding modes 
RZ, RI, and RNU. 

Injection based rounding [9] was introduced to further simplify rounding. Injection based 
rounding reduces these three rounding modes to truncation. This reduction is obtained by adding 
an injection that only depends on the rounding mode. 

A.4.2 Directed rounding of carry-save numbers 

We propose an implementation of Algorithm 1 in which intermediate results are represented by 
carry-save encoded digit strings. Since Algorithm 1 uses directed roundings for all intermediate 
computations, we need to explain how directed rounding is applied to carry-save numbers. 

Notation. Consider three binary vectors x, y, z. Let FA denote 3 : 2-addition, namely, a line of 
full-adders fed by x, y, z outputs two binary vectors s,c that satisfy \x\ + |y| + \z\ = \s\ 4- |c|. 
We regard a carry-save encoded digit string also as two binary vectors. Truncating a carry-save 
encoded digit string a[i : j] from position q simply means chopping off the tail and leaving only 
a[i : q]. We denote truncation of a from position q by \p\ q . 

We use injections to define rounding up of a carry-save digit string. The injection creates a 
carry into position q if the tail o[q -f 1 : t] is not all zeros. Formally: 

Definition 1 Rounding up of a carry-save digit string <r[0 : t] at position q is defined by 



RU q (o[0 : t]) ± [FA(a[0 : t], INJ^, q))] qi 
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where INJ^fo q) = 2~ q + l - 2"'. 

The following lemma provides bounds on the absolute error introduced by truncating and 
rounding up of a carry-save encoded digit string. 

Lemma 8 1. \a[0 : t]\ - ||*[0 : t]\ q \ e [0,2-« +1 ). 

2. \RU q (v[0:t])\-\*[0:t]\€[0 t 2-*+ l ). 

A.4.3 Dewpoint Rounding 
Dewpoints 

The concept of dewpoint rounding is inspired by the following concept of a dewpoint in meteo- 
rology: the dewpoint is the threshold temperature that separates between the events of rain and no 
rain. 

If the quotient's estimate is close enough to the accurate un-rounded quotient, then rounding 
only needs to select between two values. We denote by RC\ < RC 2 the two candidate rounded 
results. Note that RC 2 = RCi + 2~ p+1 . For each rounding candidate RC jy we denote by Ij the 
rounding interval that comprises the pre-image of RCj with respect to rounding. Namely, Ij is the 
set of numbers that are rounded to RCj. The definition of RZ, RNU, and RI rounding implies that 
since RC\ and RC2 are successive representable numbers, the rounding intervals i\ and I 2 share 
an endpoint. This common endpoint is called the dewpoint. 

To guarantee only two rounding candidates the following assumption on N' k must hold. 

Assumption 9 Prior to the computation of the dewpoint, the quotient's estimate N' k satisfies: 0 < 
p(N' k ) < 2-r. 

Dewpoint Computation 

Given a carry-save encoding a[0 : t] of N' k> we wish to compute the dewpoint and the rounding 
candidates RC\ and RC<i. A uniform computation method based on injections is used for all 
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rounding modes. 

We use some notation related to the target precision p. 

Definition 2 We refer to multiples o/2~ p+1 as representable significands and to odd multiples of 
2~ p as mid-points. 

Since 0 < p(N' k ) < 2~ p , it follows that ^ 6 [N' k , N' k 4- 2"^). A carry-save encoded digit 
string a[0 : t] is used to represent N' k . Let tail\p — 1 : t] denote the binary string that satisfies: 

\tail\p~ l:t]\ = \o\p : t]\. 

Truncation of a[0 : t] starting at position p - 1 yields: |cr[0 : p - 1]\ € (N' k - 2 • 2^ p+1 , N' k ]. By 
adding in tail\p — 1] we obtain a tighter estimate: 

MO : j> - 1]| + \tail\p - 1]| e (N k - 2~ p+1 , N' k ). 

Let RAW - \a[0 : p - 1]| + \tail\p - 1]|. We conclude that 

i^rf € [iMW, + 2 • 2" p+1 ). (A.3) 

\**\ 

Observe that RAW and RAW 4- 2~ p+1 are the only representable significands in this interval. 
This does not guarantee yet that these are the only two rounding candidates. In RI and RNU the 
interval [iL4W, RAW + 2 • 2~ p+1 ) intersects three rounding intervals. We would like this interval 
to contain at most two rounding intervals. 

Rounding mode RZ. In the RZ rounding mode the interval [RAW, RAW + 2 • 2 _p+1 ) is the 
union of exactly two rounding intervals. Hence, the dewpoint in this case is RAW + 2~ p+1 and the 
rounding candidates are simply RC\ = RAW and RC 2 = RAW + 2~ p+1 . This case is depicted 
in the top part of Fig. A.2. 
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Other rounding modes. We reduce the rounding modes RI and RNU to RZ by using injec- 
tions [1 1] as follows. 

Let I N J rn d{mode) denote an injection constant that depends only on the rounding mode. We 
define I N J rnd (mode) as follows. 



INJ rn( i{mode) := < 



0 for mode = RZ 

2~ p for mode = RNU 

2-p+i - ulp for mode = RI. 



The addition of I N J rnd (mode) reduces RI and RNU to RZ in the following sense [11]: 



RNU(x) = RZ(x -h INJ rnd (RNU)) 
RI(x) = RZ(x + INJ rnd (RI)). 

We would like to compute RC U RC 2 and the dewpoint in the cases of RNU and RI by adding 
INJ rn d(mode). One should observe that I N J rnd (mode) + \a\p:t]\ may be greater than or equal 
to 2 • 2~ p+1 . This implies that an extra bit should be used for representing the carry part of the tail. 
Namely, Let tail\p — 2 : t] denote the binary string that satisfies: 



\tail\p - 2 : t}\ = \a\p : t]\ + I N J rnd (mode) . 

Define RAW as follows: 



RAW = \a[0 : p — 1]| + \tail\p -2 :p- 1]|. 
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Claim 10 For every mode G {RZ,RNU,RI} 

RAW - INJ rnd (mode) < 7^ < RAW - INJ rnd (mode) + 2" p + 2 . 

\B\ 

Proof: Recall that \a[0 : t]\ < ^ < \a[0 : t}\ + 2"P +1 . By the definition of tail\p - 2 : t], it 
follows that |a[0 : p - 1]| + INJ rnd {mode) - 2~p +1 < jRAW* < |a[0 : p - 1]| + INJ^mode). 
The claim now follows. □ 
The proof of the following claim can be derived from Fig. A.2. 

Claim 11 For every rounding mode, RCi = RAW and RC 2 = RAW + 2 _p+1 . 

As illustrated in Fig. A.2 the difference between RC\ and the dewpoint depends only on the 
rounding mode. We introduce a second injection constant INJ dew (mode) defined as follows: 



INJ dew (mode) := < 



2-P+ 1 for mode = RZ 
2" p for mode = RNU 
0 for mode = RI, 



From Fig. A.2 it can be seen that RC\ — dewpoint = INJ^^mode), for every rounding mode. 
Hence the dewpoint is simply computed by adding RAW + I N J^imode) . 

Dewpoint back multiplication 

The rounded result Q is determined as follows: 



Q = { 



RCI if ^ < dewpoint 

RC2 if ^ > dewpoint 

RCI if ^ = dewpoint and mode = RI 

RC2 if ^ = dewpoint and mode = RZ 



(A.4) 



The case ^ = dewpoint and mode = RNU is missing since by Observation 7 it cannot occur. 



35 



The comparison <, =, > between ^ and dewpoint is done by computing the sign of 
dewpoint • \B\ — \A'\ and testing if it is zero. This task is often referred to as back multiplica- 
tion. 

A.5 Hardware Optimizations 

A.5.1 Redundant Feedback & Partial Compression 

Our FP-DIV micro-architecture is based on a Booth Radix-8 multiplier. Such multipliers are pop- 
ular in current microprocessor designs that use short clock periods. Booth Radix-8 multipliers 
are usually implemented using a 3-stage pipeline (i.e. prepare 3 A, addition tree, and final 2:1- 
addition). Goldschmidt's algorithm performs only 2 multiplications per iteration which creates 
un-utilized cycles (i.e. "bubbles"). The Athlon (in hardware) and Itanium (in software) proces- 
sors allow other multiplications to be executed during the bubbles that occur during a division 
operation. 

Instead we follow Daumas & Matula [7] and allow redundant operands so that the effective 
latency is reduced to 2 cycles. This is obtained by feeding back results represented in carry- 
save format. This design is not symmetric in the sense that the multiplier and the multiplicand 
are processed differently. The multiplier is partially compressed before being fed to the Booth 
encoder. Non-redundant representations of the multiplicand and its multiple by 3 are computed by 
non-redundant adders. 

A.5.2 Implementation of injection-based Rounding 

Injection values are fed as constants to muxes that select the injection according to the cycle and 
the rounding mode. The use of injection based rounding reduces directed rounding to truncation. 
Extra rows in the adder tree are used for adding in injections, feedback from a upper part of a long 
multiplication, or the complement of A' (during back-multiplication). 
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A.5.3 Normalization by checking (A > B) 

Suppose that the operands satisfy \A\, \B\ E [1,2). By comparing them, and shifting A to the left 
by one position if \A\ < \B\, we can normalize the quotient. Namely, we can guarantee that the 
quotient is in the binade [1,2). This normalization also slightly improves the bound on the error. 

A.5.4 Implementation of Saturated Rounding 

The definition of saturated rounding from Assumption 2 involves modifications only for the com- 
putation of F{, for i > 0. 

We propose the following implementation. During each multiplication by F{ that takes place 
in stage 2 of the pipeline, check if F' < 1. If F{ < 1, then simply de-activate the clock enable 
signal of the register that is supposed to store the result of the multiplication. The effect of ignoring 
the multiplication means that the multiplicand is not changed, i.e., precise multiplication by 1 is 
obtained. 

Note that testing if F[ < 1 can be done by testing whether the truncation of D\ is greater than 
or equal to 1. Recall that D\ is represented in carry-save format. To subtract 1 from the truncation 
of D'i we simply use a 3 : 2-adder (PPM-add) to obtain a borrow-save number. The borrow-save 
number is fed to a subtracter, and the decision is made according to the value of the sign bit. 

A.5.5 Half-size multiplication in the last iteration 

The precision in the last iteration must be at least p to obtain a sufficiently close approximation of 
the quotient. This presumably implies that the length of the second operand (i.e. the multiplier) 
should be roughly p. Based on the error analysis and, in particular, based on Corollary 6, we are 
able to implement the last iteration with a single pass through a multiplier, the length of its second 
operand being roughly p/2. 

Consider double precision. The last iteration requires computing the product N' 2 — N[ - F[. 
Recall that the length of F{ is roughly p = 53, otherwise the relative error f[ would be too large. 
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How can we avoid having to multiply N[ by the full length of F[l 
Rewrite as follows: 

N' 2 = JVJ.(1 + (JS*-1)) (A.5) 

= Ni + N[-(F[-1). (A.6) 

By Corollary 6, (F[ — 1) € [0,5i]. Hence, we may ignore bit positions to the left of position 
j = [log 2 ^J- We conclude that N£ = N[ + N[ • F[[j : t] 9 where t is roughly p. 

The addition of N[ to the product N[ • : t] can be done by using one or two of the extra 
rows in the adder tree used also for feeding injections. Note that RZ is used for the computation of 

hence these extra rows are free during this computation. 

A.5.6 Optimized implementation of Dewpoint Rounding and Back multipli- 
cation 

Recall that selecting the correct rounding candidate is based on the sign-bit (positive, negative, 
zero) of dewpoint • |J5| — \A'\. 

For the consideration of negative numbers and signs, we need to consider two's complement 
representations. In particular we consider a representation called two's complement carry-save 
representation [7]. In this representation the most significand position has a negative weight. Our 
goal in the back multiplication is to compute the sign of a = B • dewpoint — A\ we also want to 
know if a is precisely zero to determine the condition for rounding mode RI. Let z-sign(a) be two 
bits, the first signifies if a is negative, and the second bit signifies if a is zero. Hence our goal is to 
compute z-sign(a). 

We employ the multiplier for this purpose using two passes. In the first pass we compute 



a H [0 : 81] <- B[0 : 52] • dewpoint[0 : 29] + (4 - A')[-l : 52]. 
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In the second pass we compute 

a L [29 : 105] <- B[0 : 52] • dewpoint[30 : 105] + a H [30 : 81]. 

Hence, a = — a H [0] + ^[1 : 29] + a L [29 : 105]. Observe that a is represented by a two's comple- 
ment carry-save number and that with the help of a half-adder line the carry-save representations 
of a H [l : 29] and a L [29 : 83] could easily be concatenated to a carry save representation of a. It 
seems that our computation should also deal with digits in positions [—2 : —1], instead of consid- 
ering position [0] to be the most significand position. We explain now why, in fact, it suffices to 
deal only with digits in positions [51 : 105]. 

Consider a binary representation a_ 2 , a_i, . . . , aios of a, namely: 

105 

a = -a_ 2 • 2 2 + 53 a i ' 2_i - 

i=-l 

Since -2~ 52 < a < 2~ 52 , it follows that -a_ 2 * 2 2 + Z)i=-i a i * 2_i e {- 2 " 51 > °}- Hence, the bits 
in positions [—2, 51] are either all zeros or all ones. We conclude that 

105 

a = 0 iff \/ ^ = 0, (A.7) 

i=52 

a<0 iff a 51 = l. (A.8) 

A.6 Implementation & Evaluation 

This section describes the microarchitectures of the proposed division implementations. We begin 
with a discussion of the three basic stages that are involved in the computation, followed by a dis- 
cussion of the scheduling and the required operand widths for IEEE single precision and double 
precision operands. ^From the basic three stage microarchitecture we then derive several installa- 
tions, that allow to trade-off hardware cost and throughput. 
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A.6.1 Basic Microarchitecture 

A block diagram with the three stages of our basic microarchitecture is depicted in Fig. A.4. It 
arranges the blocks and techniques that were introduced in the previous sections around the core 
of a Radix-8 recoded Multiplier implementation. In the multiplier previously computed products 
can be fed back to the multiplier either as a multiplicand or as a multiplier. Latency is reduced 
to by allowing the feed back product to be a carry-save encoded digit string. The multiplier also 
supports a multiply-add operation (i.e. A x B 4- C). The addend (i.e. the number to be added to 
the product) can be input as a carry-save number or as a binary number, so 2 rows of the adder tree 
are allocated for the addend. 

The multiplication circuitry is divided in to 3 pipeline stages. Below we describe some details 
of the stages: 

1. In the first stage, the multiplicand and multiplier are prepared for the addition of the partial 
products in the second stage. The multiplier is recoded in Booth Radix-8 digits and the 
partial products are generated. Following [7], the recoding can accept either a binary string 
or a carry-save encoded digit string. For the last iteration and IEEE rounding the multiplier 
is pre-processed by a fixed shift and a redundant dewpoint computation. 

The multiplicand is processed as follows. The 3x multiple of the multiplicand is computed 
using an adder. A feed back product, encoded as a carry-save digit string, can be used as a 
multiplicand as follows. The computation of the 3x multiple is preceded by a 4 : 2-adder 
that computes a carry-save encoding of the 3x product. This carry-save encoded digit string 
is compressed to a binary number by the adder. Meanwhile, the binary representation of the 
multiplicand is computed by the adder in the third pipeline stage. This saves an adder in the 
first pipeline stage (observe that the adder in the third pipeline is indeed free). 

2. In the second stage, the partial products are compressed by a 'half-sized' adder tree. In 
addition to the partial products, there are two rows that are dedicated for (a) a carry-save 
number that is fed back from the previous multiplication, or (b) an injection (The injection 
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is used mainly to reduce a ceiling rounding to a floor rounding.). The partial products and 
the 2 extra rows are added to produce a carry-save representation of their sum. Before that 
the inputs for the two extra row are generated or selected. 

3. In the third stage, an adder compresses the carry-save number into a binary number. This 
adder can also be forced just to compute an increment. This operation mode is required to 
compute RC2 from RC\ in the last cycle of the division for rounding. Additional circuitry is 
installed to check for the condition F- < 1 from the carry-save representation of F(. and to 
check for the condition dewpoint B — A < 1 for the rounding decision from the carry-save 
representation of this expression. Few additional logic controls the selection between RC2 
from Rd based on the evaluation of the condition dewpoint • B — A < 1. 

A.6.2 Scheduling of the Stages 

In this section we present the steps of the division algorithm in detail. The division algorithm 
for single precision is listed in Table A.2 implementing one iteration of Goldschmidt's algorithm 
and the division algorithm for double precision is listed in Table A. 3 implementing two iterations 
of Goldschmidt's algorithm. The three columns of the tables correspond to the use of the three 
pipeline stages of the basic microarchitecture (although there might be more than one instance of a 
particular stage of the microarchitecture for the division implementation). As usual, we ignore the 
issues of computing the exponent and sign bit of the quotient, as these are rather straightforward. 

A.6.3 Precisions of the Operands 

There are several parameters of the implementation that need to be determined, such as the operand 
widths of the multiplier and the initial approximation. We are targeting the implementaiton for 
single precision (with p = 24) and double precision (with p = 53). Because the requirements for 
single precision are weaker than the requirements for double precision, we are only focusing on 
double precision in the following. 
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For the initial reciprocal approximation there are numerous implementation choices and one 
could even trade off initial approximation quality against internal operand precision requirements 
as it is demonstrated in Appendix C. We just assume to use the initial approximation implemen- 
tation as mentioned in [28]. This gives us |e0| < 2~ 1375 . We assume that the initial reciprocal 
estimate has a representation using less than 16 bits. 

To fulfill the requirements for dewpoint rounding according to Assumption 9, we require that 
0 < p(N' 2 ) < 2" 53 . 

With this background we are ready to determine the operand widths for the rectangular multi- 
plier implementation in our microarchitecture. Note that the multicand input A is used to feed in 
the values of N- and D[ and that the multiplier input B is used to feed in the values of F/. 

Operand Width for the Multiplier B In double precision the requirement for the widest 
operands of the multiplier B are imposed by the representation of the operands Fq and (1 — F{). 
Because operand Fq is used in carry-save representation, it requires a representation in bit positions 
[0 : \-log 2 (fo)] 4- 1]. For the operand (1 - F[) we know from Corollary 6 that (1 - F{) e [0, <Ji). 
(The case that F{ does not need to be considered because of our implementation as described in 
section A.5.4). Because (1 — F{) is available in binary (see schedule in table A.2), it requires 

representation in bit positions [L—/oy 2 (^i)J : r~'°52(/i)l]- 

For a given target approximation there are several different combinations of operand widths 
for the two operand representations that meet the requirements. Our objective in the search for 
the optimized combimation depends on the organization of our division implementation and there 
is a different objective depending on whether the two operands are succesively fed into the same 
multiplier or into two different multipliers that allow for different bit widths: 

Optimization objective 1 (minimum max bit width): 



minimumf^fote ope rand width pairs (^ox (\-log 2 (fo)] + 2, (|"-Zop 2 (/i)l - L-^2(<*i)J)) 
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Optimization objective 2 (minimum accumulated bit positions): 

\-lo92U0)] + 2 + 1 + (Hc«i(/i)l - [-log 2 (S 1 )\) 

We are only sketching the optimization for objective 1 here. Figure A.3 depicts the feasible bit 
width combinations and highlights the optimized choice regarding optimization objective 1. In the 
optimized setting the multiplier is required to be represented with 30 bits. The feasible pairs have 
been computed with the help of Theorem 4 and the target of 0 < p(A^) < 2~ 53 while assuming 
that 7r 2 only has a minor contribution. 

Operand Width for the Multiplicand A For the selected operand bit width for the multplicand, 
we need to additionally consider the requirements for 7r 2 with this setting. With Theorem 4 we 
get the requirement 7r 2 < 2~ 58,53 which implies n 0 < 2~ 60 85 making it necessary to have the 
multiplicand represented by 62 bits. 

In this way we determine our rectangular multiplier to have the dimensions 62 x 30 bits. 

A.6.4 Other Implementation Options & Evaluation 

In the previous sections the basic architecture of our division implementation has been described, 
also how the computation steps need to be scheduled to target IEEE single precision and double 
precision and how to derive the bit widths for the multiplications. 

Based on this general organization, we are proposing several different variations for the real- 
ization, that differ from each other by how many hardware components are utilized. 

Proposed Design (1/2) With a single instance of the microarchitecture from figure A.4 imple- 
menting a three stage pipeline, one can easily check from the scheduling in tables A.2 and A.3 that 
the latency for single precision is 9 clock cycles (including 1 cycle for reciprocal approximation) 
and that the latency for double precision is 11 clock cycles. This latency will also be required 
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by the other realizations that we suggest for the implementation. With one 'half-sized' rectangu- 
lar 62 x 30 bit multiplier one can schedule a new independent division every 6 cycles in single 
precision and every 8 cycles in double precision. These values of the latencies and restart times 
are listed in table A. 1 to be compared with the corresponding values from the implementations in 
current microprocessor designs. Note that for both, the latency and the startup time, our proposed 
design (1/2) shows a significant improvement, while it allows for very fast pipeline stages. Be- 
cause we only use a 'half-sized' 62 x 30 bit multiplier for this realization, the cost will only be 
roughly half of the cost of a regular multiplicative division implementation. This justifies the name 
of Proposed Design (1/2). 

Proposed Design (1), (2), and (3) By utilizing more than just one instantiation of the microar- 
chitecture from A.4, it is possible to increase the throughput of the division implementation while 
increasing the cost of the implementation. The propsed design (1) combines two 62 x 30 bit in- 
stances of the three stage pipeline to improve the restart times for single and double precision to 
3 and 4 cycles respectively. The propsed design (2) combines three 62 x 30 bit instances and one 
62 x 15 bit instance of the microarchitecture, and the propsed design (2) combines five 62 x 30 bit 
instances and two 62 x 15 bit instance of the microarchitecture. A rough cost estimation classifies 
these realizations to require about 1 x , 2 x , and 3 x the hardware cost of a conventional multiplica- 
tive division implementation, which is the reason why they are named by Proposed Design (1), 
Proposed Design (2) and Proposed Design (3). The corresponding latencies and restart times for 
the our designs are listed in Table A.l. In this way we provide a selection of options to allow for 
trade-offs between hardware cost and throughput of division implementations. The options range 
from roughly one half up to three times the cost of a regular multiplicative division implementa- 
tion, and the throughput for double precision ranges from one division every 8 cycles up to the 
fully pipelined implementation that allows to start a new division every cycle. 
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Figure A.2: Computing RCi, RC 2 and the dewpoint for the three rounding modes RZ, RNU and 
RI. The injection constant INJm d shifts RAW so that (i) the interval [RAW - INJ rnd , RAW - 
IN Jmd 4- 2~ p+2 ) contains the exact quotient, and (ii) this interval contains exactly two rounding 
intervals. 
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Figure A.6: Feasible (n, /) pairs using Setting IV for double precision are located in the shaded 
region above the curve. The closed curves depict pairs that lead to designs with equal costs. The 
central point depicts a feasible pair that leads to the cheapest design. 
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Figure A.7: Top: Feasible (n, /) pairs using setting IV for 68-bit precision (when e 0 = 2 13 51 ). 
Bottom: Cost of design as a function of /. 
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Appendix B 

Pipelined Multiplicative Division with IEEE 
Rounding 

abstract 

We propose optimized pipelined implementations for Goldschmidt's division algorithm with IEEE 
rounding based on Booth radix-8 multiplication. The considered optimizations for the quotient 
approximation are based on a careful general analysis of tight error bounds for the implementa- 
tion and are accompanied by the utilization of redundant representations, partial compressions, 
injection-based rounding and rectangular multipliers for the internal computations. For the imple- 
mentation of IEEE rounding, we introduce the concept of dewpoint rounding, that allows efficient 
implementation and reduced requirements for the quotient approximation. 

On this basis we propose the implementation of different versions of Goldschmidt's division 
algorithm with different pipeline depths. None of these implementations requires a full-sized mul- 
tiplier at any stage of the computations. In this way we reduce latency, cost, and enable increased 
throughput at a reasonable cost. We suggest a full range of pipelining depths: On one extreme is 
a 3-stage pipeline with a restart time that simply equals the latency minus the number of pipeline 
stages. On the other extreme is a fully pipelined design. 
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B.l Introduction and Summary 

Floating point multiplication and division in contemporary microprocessors is implemented by 
dedicated hardware. Moreover, designs based on multiplicative division contain a full precision 
multiplier, a rather costly circuit even in the case that it is shared for FP multiplication. In the top 
part of Table B.l, we list the latencies (i.e., number of cycles required to complete an instruction) 
and restart times (i.e., number of cycles that elapse till a functional module can engage in a new 
independent computation) of floating point division modules in several microprocessors. One can 
readily see that floating point division is rather slow compared to addition and multiplication. 

In this paper, we present implementations of Goldschmidt's floating point division algo- 
rithm [16] with significantly reduced latency and cost. We consider two settings for floating point 
division with a precision of p bits. These settings are characterized by the precision of the initial 
approximation (e.g., lookup table) and the size of the multiplier. 

• In the first setting, the precision of the initial approximation is roughly p/2 bits, and the 
multiplier is "full-sized" (i.e. the dimensions of the multiplier are roughly p x p). In the first 
setting, the latency is 9 clock cycles. 

• In the second setting, the precision of the initial approximation is roughly p/4 bits, and the 
multiplier is "half-sized" (i.e. the dimensions of the multiplier are roughly p x p/2). In the 
second setting, the latency is 11 clock cycles. 

We propose implementations for both single precision (i.e., p s = 24) and double precision 
(i.e., p d = 53). The precision of the initial approximation is roughly p d /4 « p*/2. The multiplier 
dimensions are roughly p d xp d /2& 2p s x p 3 . This implementation corresponds to the first setting 
with respect to single precision, and to the second setting with respect to double precision. 

The critical path in the designs we propose is a 62 x 30-bit rectangular multiplier that resides 
in the second pipeline stage. In fact, the multiplier is a radix-8 Booth multiplier, and the Booth 
digits of the multiplier are precomputed in the first pipeline stage. Hence, the reduction in latency 
is twofold: (i) fewer clock cycles and (ii) shorter clock periods due to a shorter critical path. 
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We propose four designs with different restart-times (but with the same latency and critical 
path). On one extreme, we suggest a three-stage pipeline with a restart-time that simply equals the 
latency minus the number of pipeline stages. On the other extreme, we suggest a fully pipelined 
design (i.e., restart-time of one cycle). These designs suggest a reasonable tradeoff between cost 
and throughput that allows for a full range of restart-times. In the lower part of Table B.l, we list 
the latencies and restart-times of the four designs we propose. 

Both extremes are quite practical. For example, in IEEE double precision, the fully pipelined 
design contains five rectangular half-sized multipliers (i.e., size 62 x 30) and two quarter-sized 
multipliers (i.e., size 62 x 15). The number of partial products associated with these multipliers is 
roughly 3.6 times the number of partial products associated with a full-sized multiplier (i.e. 56 x 
56). On the other hand, the three-stage pipeline design contains only one half-sized multiplier. 

Related work. Goldschmidt [16] developed a floating point division algorithm in which each 
iteration requires two independent multiplications. The algorithm maintains three numbers iV* 
(which converges to the quotient), Di (which converges two 1), and Fi = 2 — Z?< (called the 
scaling factor). The algorithm converges quadratically. Loosely speaking, this means that the 
number of bits needed to represent Fi doubles in each iteration. This leads to the motivation of a 
fiilly pipelined design using a sequence of multipliers (one or two per iteration), where the second 
operand's length doubles in this sequence. 

Probably the main hurdle in implementing Goldschmidt's algorithm is analyzing the errors 
due to imprecise intermediate computations. InAppendix C a general error analysis of a version 
of Goldschmidt's algorithm is presented. The advantage of the error analysis is demonstrated by 
analyzing a micro-architecture similar to that in [28] and showing that a smaller multiplier could be 
used. An error analysis similar to that in Appendix C enables one to realize the idea of a sequence 
of multipliers (one per iteration), where only the last multiplier is "full-sized". We further carry an 
idea of Goldschmidt utilizing the fact that the scaling factor tends to 1 so that we can even make 
the last multiplier "half-sized". 
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Techniques* The first step in this paper is an extension of the error analysis in Appendix C (see 
Sec. B.3). The main reason for extending the analysis is the usage of redundant representations 
(e.g. carry-save) for intermediate results. Redundant representation incurs two problems: (i) an 
increase of relative errors due to truncation (whichAppendix C easily handles), and (ii) inability to 
determine prior to compression if an intermediate result is less than (or greater than) 1. 

To reduce the drifting of intermediate results, we suggest to guarantee that the scaling factor 
Fi (which tends to 1) never falls below 1. For this purpose we introduce an assumption which 
we call saturated rounding of F { . Saturated rounding (and the fact that the scaling factor tends to 
1) enables us to use half the precision in the last iteration. This means that in the last iteration, 
instead of using a full-sized multiplier (or two passes through a half-sized multiplier), we only 
require a single pass through a half-sized multiplier. The reason that this methods works is that, 
as the scaling factor tends to one from above, there will be a run of zeros to the right of the binary 
point in the binary representation of the scaling factor. We employ our error analysis to prove an 
upper bound on the scaling factor in the last iteration. This enables us to prove a lower bound on 
the length of the run of zeros in the binary representation of the scaling factor We note that a 
similar observation (but without the one-sided convergence) was made by Goldschmidt [16]. 

Another contribution of this paper is a new rounding procedure for IEEE floating point division 
(see Sec. B.4). We refer to this rounding procedure as dewpoint rounding. The procedure relies on 
an error range of the quotient that allows for only two candidates for the final IEEE rounded result. 
We associate with each candidate r a rounding interval I r . The rounding interval I T is simply the 
set of numbers that are rounded to r. The number that separates the two rounding intervals is 
called the dewpoint. The rounding decision is obtained by comparing the dewpoint and the exact 
quotient by applying back-multiplication. We present a unified dewpoint rounding procedure for 
all rounding modes and avoid rounding tables. Previous approaches for IEEE rounding work as 
follows: (i) Compute a rounding representative of the exact quotient (i.e. a number that belongs to 
the same rounding interval that the exact quotient belongs to) (ii) Round the representative. 

We employ additional techniques, among them: (i) A Booth recoded multiplier that can be fed 
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by either non-redundant binary operands or by redundant carry-save operands [7]. This technique 
when applied to a Booth radix-8 multiplier enables reducing the feedback latency to two cycles. 

(ii) Injection based rounding is used to implement directed rounding of intermediate results [9]. 

(iii) Back-multiplication (i.e comparison between the dewpoint and the exact quotient) requires a 
full precision multiplication. In the design that uses a half-sized multiplier, two passes through 
the multiplier are required. To avoid wasting a cycle in preparing the exact dewpoint, we use an 
approximation of the dewpoint in the first pass (this reduces the latency because, otherwise, we 
would insert a bubble during this cycle into the pipeline). In the second pass, a correction is used, 
so that the error caused by the first pass is fully corrected. This method of using an approximate 
dewpoint heavily relies on the error analysis. 

Organization. In Section B.2, notation is introduced, division is defined, and the version of 
Goldschmidt's algorithm from Appendix C is reviewed. An error analysis of this division algo- 
rithm appears in Section B.3. This analysis builds on the analysis that appears in Appendix C. In 
Section B.4, dewpoint rounding is presented. In Section D.6, we present an implementation of the 
division algorithm that is based on a full-sized multiplier. For the sake of concreteness, the presen- 
tation is for single precision and the implementation uses a double-precision micro-architecture. 
This means that the half-sized multiplier for double precision is a full-sized multiplier for single 
precision. In Section B.6, we present an implementation of the devision algorithm that is based 
on a half-sized multiplier. For the sake of concreteness, an implementation for double precision is 
presented. In Section B.7, we list various techniques employed to save cost and delay in a hard- 
ware implementation. In Section B.8, we consider three pipelining options ranging from a fully 
pipelined design to a design that reduces the restart-time by half. 
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Table B.l: Latencies(restart-times) of floating-point operations in current commercial micropro- 
cessors compared to the proposed division implementations. 

B.2 Preliminaries 

Notation Let XiX i+ \ • • • xj G {0, 1}* denote a binary string. We often denote this string also by 
x[i : j]. We also sometimes refer to x» as x[i}. Since we deal with fractions (mostly in the binade 
[1,2)), the weight associated with the bit Xi is 2~\ Namely, a fraction is represented by the binary 
digit string x 0 .xix 2 • • • x p - X . 

Let a = cricr i+ i • • • aj G {0,1, 2}* denote a carry-save encoded digit string (i.e. a* G {0,1, 2}). 
A carry-save encoded digit string x[i : j] is represented by two binary vectors xc[i : j] and xs[i : j]. 
Each carry-save digit x[£] satisfies x[£] = xc[£] +xs[£]. We refer to xc (resp. zs) as the carry-string 
(resp. sum-string) of x. 

Given a binary string j] and a carry-save encoded digit string a[i : j], we denote the values 
represented by these string by \x[i : j] \ and |cr[s : j]\, respectively. 
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Division Operation We consider the following task which captures the main difficulty in IEEE 
floating-point division. The dividend and the divisor are represented by p-bit binary strings A[0 : 
p - 1] and B[0 : p - 1], where \A\, \B\ e [1, 2). Our goal is to compute Q[0 : p - 1] that is the 
binary encoding of the rounded value of the normalized quotient. More formally, (i) Let A 1 = A 
if \A\ > \B\, and A' = leftshift(A) if \A\ < \B\. This way, the quotient \A'\/\B\ € [1, 2). (ii) let 
q = (iii) round q (according to the IEEE standard). The binary string Q[0 : p — 1] should 

satisfy |Q| = q. 

We consider two precisions: double precision in which p d = 53 and single precision in which 
p s = 24. Four rounding modes are defined in the IEEE standard: round towards zero round 
towards zero (RZ), round to nearest even (RNE), round towards infinity (RI) and round towards 
minus infinity(RMI). 

Goldschmidt's division algorithm InAppendix C a variation of Goldschmidt's division algo- 
rithm based on directed rounding is presented and analyzed. The algorithm is listed below as 
Algorithm 2. The multiplications involved and their dependencies are depicted in figure B.l. 

Directed roundings are used for all intermediate calculations. For example, iV t ' is obtained by 
rounding down the product N!_ x • F-_ v We denote by the relative error incurred when iV-_i • 
is rounded down. Rounding down translates to the assumption that > 0. Similarly, rounding 
down is used for computing F( and rounding up is used for computing D\. Our error analysis 
requires that: (a) the operands are in the range A\ B satisfy A e [B, 2B) and B € [1, 2), so that 
A'/B e [1,2) (b) all the relative errors incurred by directed rounding are at most 1/4, and (c) 
we require that |e 0 | + 3d 0 /2 + / 0 < 1/2. (Note that neither of these requirement imposes any 
restriction and is very easily met.) 
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D',F'-Pipeline JN'-Pipeline 
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Figure B.l: Schedule of the Iterations of Goldschmidt's division algorithm using approximate 
arithmetic and an initial approximation for l/B. The numbers in circles indicate the sequence 
of the multiplications involved. For our implementation a bound on the relative error p{N[ ) of 
iteration i appears in each iteration. 



B.3 Extension of error analysis 

In this section we present a variation of the error analysis presented in Appendix C. The error 
analysis presented in Appendix C used a simplifying assumption called strict directed rounding 
(SD-rounding). Satisfying SD-rounding is easy if intermediate results are represented in non- 
redundant binary representation. In our division algorithm redundant representation is used for 
intermediate results (i.e. carry-save and borrow-save encoded digit strings). Hence we cannot rely 
on the SD-rounding assumption. 

The analysis uses the same requirements used before, namely: (a) the operands A' and B satisfy 
A! e [B, 2B) and B 6 [1, 2), (b) all the relative errors incurred by directed rounding are at most 
1/4, and (c) |e 0 | + 3d 0 /2 + /o < 1/2. 
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Algorithm 2 Goldschmidt-Approx-Divide(A', B) - Goldschmidt's division algorithm using ap- 
proximate arithmetic 

1: Initialize: *- A', D f _ l <- B, F'_ x *- 

2: for £ = 0 to k do 



6: end for 

7: Return^') 

B.3.1 Rounding assumptions 

The strict-directed rounding assumption is defined as follows. 

Assumption 12 fSD-roundingj Directed rounding up (resp. down) r is strict ifx<l implies 
r{x) < 1 (resp. x > 1 implies r(x) > 1). 

Instead of assuming that SD-rounding holds, we use the following assumptions for iterations i > 0: 
Assumption 13 (saturated rounding; Fori > 0, 



Note that the definition of the relative error /* in the computation of F- is now ambiguous (i.e. 
suppose F( = 1 and (1 — /») • (2 — £>•) < 1). To avoid ambiguity, define the relative error // as 
follows: 



Note that if D\ < 1, then /» > // > 0. This means that saturated rounding reduces the relative 
error in the computation of F? when D\ < 1. 

The following assumption requires that multiplication by one is precise. 

Assumption 14 (precise multiplication by i; IfF[ = 1, then D' i+l = D\ and A^ +1 = N[. 



3: Nl^-il-^.NU-FU 
5: *-(!-/,)• (2-D<). 



l? = max{l, (1-/0 -(2 -£>?)}. 
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We say that the algorithm is stuck starting with iteration i if Dj — D[ and Afj = for every 
j > i. Note however, that F? = 1 does not imply that the algorithm is stuck starting with iteration 
i. The algorithm may not be stuck if D\ < 1 and eventually the relative error fj is small enough so 
that F( > 1. However, ii D[ > 1, then the algorithm is stuck starting with iteration i. 

We denote the intermediate values obtained with saturated rounding and precise multiplication 
by 1 also by N[, D\ y and F[. 

B.3.2 Analysis 

Our analysis shows that our weaker assumptions do not affect the bounds on the errors fromAp- 
pendix C. 

Theorem 15 For every i > 0, the relative error p(iV/) = A ^Jq % satisfies 

o< p(jvj) < + (B.i) 

where ix* is defined by Hi = 1 — (1 — n*) • n}=o IT^ — 0 w/iere 4 w defined by 




|e 0 |+3do/2 fori = 0 
5f_ l + /i-i otherwise 



for i > 0. 

Lemma 16 The following three statements hold: 

a) IfFU = (1 - • (2 - Z?U) ^-i € [1 - 1 + ften H{ > 1 - 
W //i^-i = (1 - /t-i) ' (2 - Z3U) then D< < 1 + d*. 

c) //i^ = 1 > (1 - • (2 - ZT^) and G [1 - 1], rten € [1 - S u 1]. 
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Proof: Let tj = - 



1. We first prove Parts a) and b): 



D\ = (1 + dj) • (1 + ti-i\ - (1 ~ • (1 ~ /.-i) 




>i <i <i 



It follows that > (1 - 6f_ x - = 1 - S u which implies Part a). It also follows that 

D\ < 1 + du which implies Part b). 

c) Since F{_ x = 1, precise multiplication by 1 implies that D[ = D[_ l9 so all we need to show 
is that D't > 1 - 8i. Since 1 > (1 - • (2 - we can write: 



Proof of Theorem 15: InAppendix C Equation (B. 1) is proven for iteration i whenever € 
[1 — 1 + 5i-i] and F(_ x = (1 — • (2 — In fact, these two conditions are used in 

Appendix C to prove that 1 — Si < Z^_x • F(_ x < 1, which implies Equation (B.l) for iteration i. 
We rely on the proof presented in Appendix C and deal here with the effect of saturated rounding. 
Let %d and ip be defined as follows: 



Dl > (1 - /<_!) • (2 - DU) ■ Dl 
>(!-/«_!)• (l-^i) 



'i-i 



> 1 -& 



and Part c) follows. 



□ 



i D = min{i >Q:D' i >\) 

i F = min{z > 0 : (1 - /•) • (2 - Z>|) < 1} 



Note that if > 1, then (1 - f { ) ■ (2 - E%) < 1, and hence, i F < i D . 
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We need to consider the 3 cases: (i.) i < %f\ ("0 i>F < i < Id\ and (iii.) i > i D . Note that case 
(ii.) only occurs \fi F ^i D , 

Case (i.): One can directly verify that 1 - S 0 < D' Q < 1 + 5 0 . Since saturated rounding does not 
take place in iteration i = 0, it follows that F£ = (1 - / 0 ) • (2 - D' 0 ). Hence the error analysis in 
Appendix C implies that Equation B.l holds for i = 1. 

For 0 < i - 1 < i F we have D^ x < 1 and F{_ x = (1 - fr-i) • (2 - D[_ x ) > 1 since 
i — 1 < min(iD> If). By applying induction and Part a) of Lemma 16, it follows that 6 [1 — 
5i-i, 1]. Hence the proof in Appendix C applies, and Equation B.l holds for i. It follows that 
Equation B.l holds for every i € [1, 

Case (ii.): By applying induction and Parts a) and c) of Lemma 16, it follows that Z^_ x 6 [1 — 
8i- u 1], for every i - 1 < i D . If FJLj = (1 - f^i) - (2 - D^), then Equation (B.l) holds for i. 
If = 1 > (1 - • (2 - !>•_!), then by Part c) of Lemma 16, D f { e [1 - 1]. By precise 
multiplication by 1 it follows that D[ = D[_ l9 and hence 1 - S { < • < 1. Hence the error 
analysis in Appendix C is applicable, and Equation (B.l) holds for i < io. 

Case (iii.) (i > Z£>): Note that the algorithm is stuck starting from iteration i D . Since Wj+i > 7r J5 
it suffices to prove that 0 < p(N[ D ) < 7i(i D +i). 

We now prove that D io € [1,1+ <k D ]. The lower bound follows from the definition of i D . 
Part b) of Lemma 16 proves that if F iD . x = (1 - f lD - X ) • (2 - D Id - X ) 9 then D iD < 1 + d iD . So 
all we need to prove is that F( D ^ = (1 - /i D -i) ■ (2 - A' D -i)« Since A D > 1 it follows that 
if A 0 -i < 1, then F[ D _ X ± 1, hence F[ D _ X = (1 - /^-O • (2 - I^). The minimality of i D 
implies that if A^-i > 1, then i D = 1. But saturated rounding is not used in iteration i = 0, so 
the upper bound on Z^ D holds. 
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JVJ 



We now expand as follows 



Thus, 



and 



<=0 * 



Since D< D > 1, it follows that 



<i-d --*,«,) n^f 



1 - n e (B.3) 

£=0 



which gives the desired upper bound on p{Nl D ). To prove the lower bound we use D\ D < 1 4- d io 
with equation (B.2) as follows: 



%D i ^ 
1 — Tit 



-i-a-^j n £s 



= 7T iD > 0. 
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This completes the proof for the last of the three cases. □ 
Corollary 17 For every i > 0, F( € [1, 1 + <*<] 

Proof: It follows from Assumption 13 (saturated rounding) that F[ > 1. We only need to deal 
with the case that F[ = (1 - /») • (2 - In the proof of Theorem 15 we showed that D\ > 1 - S i9 
for i > 0. Therefore, F( = (1 - /») • (2 - < 1 + □ 



B.4 DEEE rounding of the quotient 

In this section we describe how IEEE rounding is done in our algorithm. We present a novel 
method called dewpoint rounding. 

B.4.1 Background on IEEE rounding 

The IEEE Floating-Point Standard 754 [17] requires the implementation of all four IEEE round- 
ing modes for all basic arithmetic operations and for all precisions supported. IEEE rounding is 
supposed to be computed on the exact result of the operation (which is A'/B in our case). For 
most operations the exact result is not available or too expensive to be computed. It is therefore 
common practice to compute the IEEE rounding result by computing a rounding representative 
first and compute IEEE rounding on the representative that leads provably to the same result [10]. 

In the case of multiplicative division it is even difficult to find a representative quotient from an 
approximation of the quotient. In [18], it is proven that the length of runs of zeros (or ones) in the 
binary representation of the exact quotient is limited. This means that a rounding representative can 
be decided from an approximate quotient if the precision of the approximate quotient is roughly 
twice the target precision. We do not follow this approach since it involves one additional iteration 
for multiplicative division algorithms, adds significant delay to the implementation, and requires 
very wide intermediate operand representations for the last iteration. 
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Back-multiplication is the common method for computing the IEEE rounded quotient. In this 
method, the approximate quotient is multiplied by the divisor jB, and then this product is compared 
with the dividend A' [19, 34]. 

In the sign-magnitude representation of floating-point numbers, the four DEEE rounding modes 
can be reduced to three rounding modes RZ, RI, and RNU (round to nearest up) based on the 
numbers' sign [32]. 

Observation 18 The exact quotient A' /B cannot be a midpoint between two representable num- 
bers. Therefore RNU and RNE are equivalent rounding modes with respect to division. 

Based on the above reductions and Observation 18, we only need to consider the three rounding 
modes RZ, RI, and RNU. 

Injection based rounding [9] was introduced to further simplify rounding. Injection based 
rounding reduces these three rounding modes to truncation. This reduction is obtained by adding 
an injection that only depends on the rounding mode. 

B .4.2 Dewpoint Rounding 
Dewpoints 

The concept of dewpoint rounding is inspired by the following concept of a dewpoint in meteo- 
rology: the dewpoint is the threshold temperature that separates between the events of rain and no 
rain. 

If the quotient's estimate is close enough to the accurate un-rounded quotient, then rounding 
only needs to select between two values. We denote by RC X < RC 2 the two candidate values for 
the IEEE rounded result. Note that RC 2 = RCi + 2~ p+1 . For each rounding candidate RC jy we 
denote by Ij the rounding interval that comprises the pre-image of RCj with respect to rounding. 
Namely, Ij is the set of numbers that are rounded to RCj. The definition of RZ, RNU, and RI 
rounding implies that since RCi and RC 2 are successive representable numbers, the rounding 
intervals Ji and I 2 share an endpoint. This common endpoint is called the dewpoint. 
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To guarantee only two rounding candidates the following assumption on N k must hold. 

Assumption 19 Prior to the computation of the dewpoint, the quotient's estimate N k satisfies: 
0 < p(N' k ) < 2-P 

Computation of rounding candidates 

Given a carry-save encoding <j[0 : t] of N ki we wish to compute the dewpoint and the rounding 
candidates RC\ and RC2- A unified computation method based on injections is used for all round- 
ing modes. We begin by dealing with RZ rounding and then reduce the two other rounding modes 
RNU and RI to RZ. 

We use some notation related to the target precision p. 

Definition 3 We refer to multiples o/2~ p+1 as representable significands and to odd multiples of 
2~ p as mid-points. 

RZ rounding. Since 0 < p(N' k ) < 2~* and < 2, it follows that e [N' k , N' k + 2-P+ 1 ). 
A carry-save encoded digit string o[0 : t] is used to represent N k . Let bin(a)[0 : t] denote the 
binary string that satisfies \bin(a)[0 : t]\ = |cr[0 : t]\. (Note that this notation does not imply full 
compression of a in the implementation. We use this notation just to simplify the description of 
the rounding.) 

Truncation of bin(a)[0 : t] after position p — 1 yields: 

\bin(cj)[0 :p-l}\ e(N' k -2-* +1 ,N' k ). 

Let r = \bin(o)[0 : p — 1]|. We conclude that 

i^f €[r,r-f2.2- p+1 ). (B.4) 
1^1 

In the RZ rounding mode the interval [r,r + 2 • 2~ p+1 ) is the union of exactly two rounding 
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intervals: [r,r+2~ p+1 ) and [r + 2~ p+1 ,r + 2-2~ p+1 ). Hence, the dewpoint in this case is r + 2~ p+1 . 
The rounding candidates are simply RCi = r and RC2 = r + 2~ p+1 . 

This method fails for RI and RNU since in these rounding modes the interval [r, r + 2 • 2~ p+1 ) 
intersects three rounding intervals. 

Reduction of RI and RNU to RZ. We reduce the rounding modes RI and RNU to RZ by using 
injections [11] as follows. 

Let IN J rnd (mode) denote an injection constant that depends only on the rounding mode an on 
the precision p. We define iNJ rn d(mode) as follows. 



IN J rn d (mode) := < 



0 for mode = RZ 

2- p for mode = iWt/ (B.5) 

2-P+1-2-* for mode = i?/. 



For x E [1, 2), the addition of iNJ rn(i (mode) reduces RI and RNU to RZ in the following sense 
(c.f.[ll]): 

RNU(x) = RZ(x + \N5 rnd (RNU)) 

(B.6) 

RI(x) = RZ(x + im rnd (RI)). 

Hence it suffices to compute RZ(^ -h ^1^(77100^6)). 

We now repeat the analysis with lNJ rnd (mode) as follows. Since 0 < p(N' k ) < 2" p , it follows 
that ^ + lNJ rnd (mode) € [N' k + lNJ rnd (mode), N' k + IN J rnrf (mode) + 2" p+1 ). Let bin(a')[Q : t\ 
denote the binary string that satisfies 



\bin(c')[Q :t]\ = \a[0 : t}\ + I NJ rnd (mode). 



(B.7) 



69 



Truncation of bin(a')[0 : t] gives the following estimate: 



\Un{o')[Q : p - 1]| € (N' k + INJ rnd (mode) - 2" p+1 , iV£ + INJ rnd (mode)]. 



Let r' = |6m(<j')[0 : p - 1]|. We conclude that 



1*1 



+ lNJ rnd (mode) € [r', r' + 2 - 2" p+1 ). 



(B.8) 



Definition 4 The rounding candidates are defined as follows: 



RC X = r' 



and 



RC 2 = r' + 2" p+1 . 



The following claim summarizes the reduction. 

Claim 20 For every rounding mode mode G {RZ,RNU,RI}, the IEEE rounded result satisfies the 
following equation: 



\A'\ < (RC 2 - lNJ rnd (mode)) • 
Dewpoint back-multiplication 

The disadvantage in Claim 20 is in the complexity of computing the comparison \A'\ < (RC 2 — 
iNJ rm i(morfe)) • \B\. Specifically, in RI mode, the representation of lNJ rn d(mode) requires t bits 
of precision. 

In this section we show how back-multiplication can be done without having to deal with 
positions to the left of position \p]. 





(B.9) 



otherwise. 



v 



Note that the decision does not require division since 



+ INJ rmi (mode) < RC 2 if and only if 
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By Claim 20, the rounded quotient equals RCi if and only if 

i^i < t' + 2-<»-V - im rnd (mode). (B.10) 

Since we are interested in simplifying the computation of the dewpoint, we define the dewpoint 
displacement constant (^^(mode) as follows. 

2-< p - x > if mode =RZ 
Cdevjimode) = { i~v if mode =RNU 
0 if mode =RI. 

Note that C^ w {mode) < 2" p+1 - WJ rnd (mode) for every rounding mode. 
Definition 5 The dewpoint is defined as follows: 

dewpoint = r' + C&> w {rnode). 

Note that dewpoint = r' + 2~ (p ~ 1 ) — iNJ rnd (roode) if the rounding mode is RZ or RNU. If the 
rounding mode is RI, then dewpoint is too small by 2"*. Since 2~ l is the our resolution, it follows 
that ^ < r> + 2"^- 1 > - iNJ rnd (mode) if and only if ^ < r' + 2"^- 1 > - iNJ rnd (mode) - 2"*. 
This implies in RI, the rounded quotient equals RC\ if and only if ^ < dewpoint (as opposed to 
the strict inequality in Equation B.10). 

Back-multiplication is implemented by computing the sign (positive, zero, negative) of the 
expression 

\A'\ — \B\ • dewpoint. 

The following claim summarizes how back-multiplication is used to determine the rounded quo- 
tient. 

Claim 21 The IEEE rounded quotient Q equals RC\ or RCi. The decision is determined as 
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follows: 



RC! if(\A'\ 



\B\- dewpoint) < 0 



Q= < 



RC 2 if(\A'\ 



\B\ • dewpoint) > 0 



RCi if(\A'\ 



\B\ ■ dewpoint) = 0 and mode = RI 



RC 2 if(\A'\ 



\B\ • dewpoint) = 0 and mode = RZ 



v 



The case ^ = dewpoint and mode = RNU is omitted since by Observation 18 it cannot occur. 



B.5 An implementation with a full sized multiplier 



In this section we present an implementation of our algorithm that uses a full sized multiplier (i.e., 
for a target precision p, the multiplier dimensions are roughly p x p) and an initial approximation 
that whose precision is roughly p/2. To be concrete, we present detailed design for single precision 
(i.e. p s = 24) that has a latency of 9 cycles. 

The main motivation for this setting is for performing single precision divisions with a design 
built for double precision. We assume that core of the design for double precision division is a 
half-sized multiplier (i.e., somewhat larger than a p d x (pd/2) multiplier, where p d = 53). Such a 
multiplier is larger than a 2p s x p s multiplier. This implies that, when performing single precision 
division, the multiplier is not only a full-sized multiplier, but a double-width multiplier. We can 
employ a double-width multiplier to reduce to latency from 9 cycles to 8 cycles (see item 1 in 
Section D.6.4). 

B.5.1 Basic microarchitecture 

A block diagram with the three stages of our basic microarchitecture is depicted in Fig. B.2. The 
core of this microarchitecture is a radix-8 Booth recoded Multiplier. The microarchitecture allows 
for feeding of previously computed products back to the multiplier as either operands. Latency is 
reduced by allowing the feedback to be in redundant representation. The multiplier also supports 
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a multiply-add operation (i.e. i'xB + C). In this section the addend (i.e. the number to be added 
to the product) is a binary number, so only one row in the adder tree is allocated for the addend. 

The multiplication circuitry is divided into 4 pipeline stages. Below we describe some details 
of the stages: 

0. In pipeline stage 0, an estimate for the reciprocal 1/B is computed and the values of A and 
B are compared to determine the pre-normalized value A \ Additional hardware is included 
for scheduling and making available the operand values A', B and approx(l /B) to the inputs 
of the first and the second stage. 

1. In the first stage, the two operands of the multiplier (i.e., the multiplicand and multiplier) 
are prepared for the addition of the partial products in the second stage. The second operand 
(i.e., the multiplier) is recoded in Booth radix-8 digits and the partial products are generated. 
Following [35], the recoding can accept either a binary string or a carry-save encoded digit 
string. 

The first operand (i.e., the multiplicand) is processed as follows. The 3x multiple of the 
multiplicand is computed using an adder. A feedback product, encoded as a carry-save digit 
string, can be used as a multiplicand as follows. The computation of the 3 x multiple is 
preceded by a 4:2-adder that computes a carry-save encoding of the 3x multiple. This carry- 
save encoded digit string is compressed to a binary number by the adder. Meanwhile, the 
binary representation of the multiplicand is computed by the adder in the third pipeline stage. 
This saves an adder in the first pipeline stage (we need to insure that the adder in the third 
pipeline stage is indeed available). 

2. In the second stage, the partial products are compressed by an adder tree. In addition to 
the partial products, there is a row that is dedicated for (i) a carry-save round-up rounding 
injection (see Section D.5.2), (ii) an IEEE rounding injection (see Equation D.6), or (iii) a 
two's complement number (see Section D.5.3). 
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3. The third stage contains an adder and a comparator. The adder has the following tasks: 

(i) compressing the multiplicands from carry-save representation to binary representation 

(ii) computing RC\ and RC^ 

The comparator is used for computing the sign of the back-multiplication and checking for 
F[ < 1 from the carry-save representation of D[. We cannot use the adder as a comparator 
due to conflicts between two pipelined divisions. 

Additional circuitry is installed in the third stage, including: (i) circuitry for feeding the adder 
with 2~ 23 , (ii) circuitry for computing RCi, RCi, and (iii) circuitry for selecting between 
Rd and RC 2 . 

B.5.2 Pipeline schedule 

Table B.2 lists the schedule of operations that implements our algorithm for single precision based 
on a "full-sized" multiplier and a "half-sized" initial approximation. In the error analysis presented 
in Section D.6.3 we show that a 30 x 27 rectangular multiplier together with an error bound of 
|e 0 | < 2 13 74 guarantees correctness. We use these parameters in the detailed description below. 
We now describe the cycles one by one. 

cycle 0: In cycle 0 an initial approximation of 1/B is computed. We denote this approximate 
reciprocal by (1 — eo)/B. Since this reciprocal is recoded in the next cycle, it is possible 
to have the reciprocal represented in binary or redundant format. Moreover, in cycle 1, very 
little processing of this reciprocal is required. This means that one could actually allow the 
reciprocal approximation to take place also during part of cycle 1. 

cycle 1: Iteration 0 of the algorithm begins in this cycle. The approximate reciprocal (1 — e 0 )/B 
is assigned to FL X and is fed as the second operand of the multiplier. The approximate 
reciprocal is recoded as Booth radix-8 digit string. The first operand of the multiplier D'_ x is 
simply assigned the value B. The 3x multiple of D'^ is computed in this cycle. 
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In the third pipeline stage, A and B are compared. If A < B, then we set A' <— 2 - A 9 
otherwise A! <— A. This is how we guarantee that A' e [B, 2B). Note that since A' can be 
in the binade [2, 4), A' has a bit in position [—1], 

cycle 2: Two activities take place in this cycle. In the first pipeline stage, A'[—l : 28] is assigned 
to NL\ and the 3x multiple of N'_i is computed. In the second pipeline stage the rounded up 
product D'_\ - FL-l is stored in D' o [0 : 28] as a carry-save number. The rounding up of a result 
represented in carry-save representation is done by using a rounding up injection particular 
to carry-save representations INJCS /h/(28 } 28 4- 14). The details of this method are explained 
in Section D.5.2. 

cycle 3: The carry-save representation of D f 0 is compressed to binary in the third pipeline stage. 
In the second pipeline stage, the rounded down product • FL± is stored in N^[-l : 28]. 
Rounding down is simply achieved by truncation of the carry-save representation. Note that 
our error analysis (Theorem 15) only implies that N[ < 2 for i > 0. Hence we need position 
[-1] mN(>. 

cycle 4: In the first pipeline stage, a carry-save representation of Fq = 2 — D' 0 is computed and 
then Booth recoded. (Recall that Dq is compressed to binary in the previous cycle.) In 
addition, the 3x multiple of Nq is computed. 

In the third pipeline stage, the carry-save representation of Nq is compressed to binary. 

cycle 5: In this cycle we compute in the second stage N" in carry-save representation with the 
value: 

N" <- No * + lNJ rnd (mode). 

The carry-save representation of N" equals the value of \bin(a') \ defined in Equation B.7. 

cycle 6: In the third pipeline stage, RC\ is computed by compressing N" and truncating at position 
23. The dewpoint equals RCi + Cd^imode), However, to reduce delay, we compute the 
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dewpoint as follows. A carry-save representation of the dewpoint is obtained by adding 
N{ f [0 : 23] + Cdevimode) + carry 2 z{N'{ '[24 : 28]), where carry 2 3 (N{'[24 : 28]) is a bit that 
denotes if N"[2A : 28] > 2 -23 . The dewpoint is Booth recoded and the 3x multiple of B is 
computed. 

cycle 7: Back-multiplication takes place in the second pipeline stage. We use a representation 
called two's complement carry-save representation [7]. In this representation the most sig- 
nificant position has a negative weight. Our goal in the back multiplication is to compute the 
z-sign of a = B • dewpoint — A 1 . (The z-sign indicates equality, greater than, or less than.) 
In Section D.5.3 we show that it suffices to examine digits in positions [22 : 47] to determine 
the z-sign of a. 

cycle 8: In the third pipeline stage, (i) RC 2 is computed, (ii) the z-sign of the back-multiplication 
is determined from a[22 : 47], and (iii) the IEEE rounded quotient is selected from RC\ and 
RC 2 according to the rounding mode and the z-sign of a. 

B.5.3 Correctness 

In this section we prove the correctness of a design based on a rectangular 30 x 27-multiplier and 
an initial approximation with a relative error eo < 2~ 13 74 . To support single precision (p s = 24), 
dewpoint rounding requires that Assumption 19 holds. Namely, 0 < p(N{) < 2~ 24 . 

The analysis presented in this section is a concrete analysis for single precision. More gener- 
ally, this analysis applies to the following situation: (i) the accuracy of the initial approximation is 
roughly 2" p/2 ; (ii) the quotient is approximated by N[ (namely, the algorithm has 1 iteration of 
quadratical convergence); and (iii) the multiplier is a rectangular multiplier whose dimensions are 
not smaller than roughly p s x p s . 

In Table B.2 we assumed that initial approximation FL X = is represented by 15 bits. This 
implies that the multiplier can easily accept F'_i as a short operand. 
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Claim 22 If the relative error in initial approximation satisfies |e 0 | < 2~ 1374 , then correct double 
precision IEEE rounding is obtained with a rectangular 30 x 27-bit multiplier. 

Proof: All we need to show is that 0 < p(N[) < 2~ 24 . By Theorem 15, 0 < p(N[) < ir x + 5 X . 
We expand tti and 6i to obtain: 

1 - n 0 

7T! = 1 - (1 - m) • ■ 

/ 3 V * <""> 
<5i = ( |e 0 | + - -rf 0 ) +/o- 

To prove 7Ti 5i < 2~ 24 , we prove upper bounds on the relative error terms (i.e., n iy f iy and 
substitute these upper bounds in Equation B. 1 1 . We bound these relative error terms using absolute 
error terms defined as follows: 

Definition 6 The absolute errors of intermediate computations are defined as follows : 

nepsi = Ui • N'^ • F(_ x 
depsi ± di-D'^-FU 
fep Si = fi'V-D't). 



We now bound each of the relative error terms: 

Bound on do: According to Table B.2, D' 0 is represented by a carry-save encoded digit string with 
digits in positions [0 : 28]. This digit string holds the rounded up value of D'_ x • F!_i- This 
implies that the absolute error depso associated with the computation of D' 0 is in the range 
[0, 2 • 2~ 28 ). Since the precise product D'_ x • FL\ equals 1 - e 0 , this product is in the range 
1 ± 2" 13 - 74 . It follows that d 0 < 2" 27 /(l - 2" 13 - 74 ) < 2" 26 " 989 . 

Bound on n 0 : As in the case of D' 0 , Nq is also represented by a carry-save encoded digit string 
with digits in positions [0 : 28]. This digit string holds the rounded down value of N'_ x • FL\* 
This implies that the absolute error neps 0 associated with the computation of D' 0 is in the 
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range [0, 2 • 2" 28 ). The product A! • FL X equals § • (1 - |e 0 |). Since A' > S, it follows that 
A' • > 1 - |e 0 |, and therefore, n 0 < 2" 27 /(l - |e 0 |). 

Bound on fo: Note that Dq is compressed to binary non-redundant representation in the third 
cycle. We first show that feps 0 < 2~ 26 . The reason is that instead of computing 2 — 
D' 0 [0 : 28], we compute 2 - D' Q [0 : 26]. The absolute error is -£>o[27 : 28]. Hence, 
fepso < 2~ 26 , as required. Next, we expand 2 - D' Q = 2 - (1 -h d 0 ) • (1 - e 0 ). We conclude 
that 

2~26 

/o< 2 - (1 + db) - (1 + |co|) 

< 2~25.99989454 



Bound on m: As in the case of n 0 , nepsi < 2 27 . 

We bound from below and Fq as follows: = (1 -n 0 ) • (1 - e 0 ) • ^f. Hence, > (1 
n 0 ) • (1 - |e 0 |). Also, Fq = (1 — /o) • (2 - Dq) > (1 - fo) ■ (1 - *o). We conclude that 



2 -27 



< 2 



(l-no)-(l-|cd|)-(l-/o)-(l-ft) 

-26.999789 



We now combine these bounds to obtain: 



7*+* <2- 24 48458 <2" 24 



□ 
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B.5.4 Design alternatives 

The proposed implementation is not the only way to implement our algorithm. Certain choices 
were made in order to keep the implementation simple without increasing cost by much. Below 
we list a few alternatives. 

1. If the length of the first operand of the multiplier is wider than 29 + 30 bits, then cycles 2 
and 3 could be executed simultaneously. The reason is that the two multiplications that take 
place in these cycles have the same second operand (i.e., F!_i)- (Compression of D' Q is not 
crucial and could be skipped at the expense of increasing the error term / 0 .) This implies a 
reduction in the latency from 9 cycles to 8 cycles. 

2. The normalization of the quotient can be postponed till after N[ is computed. This modifica- 
tion has the following implications: (i) A' is not computed, and we simply use A instead, (ii) 
In this version the rounding injection cannot be added in cycle 5 since normalization changes 
the value of the injection. Hence in cycle 6 an injection (or injection correction) needs to 
be added when computing RC\. (iii) The error analysis should be changed to deal with this 
case (slightly worse bounds will be obtained). 

3. We chose to use a multiplier in which both the first operand and the product have a digit 
in position [— 1]. This can be avoided by normalizing the operand and product so that the 
leftmost digit is always non-zero. 

4. If the initial reciprocal approximation has a one-sided error, namely e 0 > 0, then we can 
also guarantee that N$ < 2. This means that the the first operand of the multiplier need not 
support position [-1]. In the second pipeline stage in cycle 3, we could compute instead: 
mulRziNLAO : 28], FijO : 14]) + NL X [-1] * FL x [0 : 14]. 
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B.6 An implementation with a half-sized multiplier 

In this section we present an implementation of our algorithm that uses a half-sized multiplier (i.e., 
roughly p x p/2 for a target precision p) and an initial approximation whose precision is roughly 
p/4. To be concrete, we present detailed design for double precision (i.e. p d = 53) that has a 
latency of 11 cycles. 

The microarchitecture is an extension of the microarchitecture introduced in Section D.6. 1. The 
main differences are in the computations related to the last iteration and the back-multiplication. 
These changes influence circuitry needed for selection, injections, computation of the rounding 
candidates, and the dewpoint. 

B.6.1 Pipeline schedule 

Table B.3 lists the schedule of operations that implements our algorithm for double precision based 
on a half-sized multiplier and a quarter-sized initial approximation. In the error analysis presented 
in Section D.6.3 we show that a 62 x 30 rectangular multiplier together with an error bound of 
|e 0 | < 2 13 74 guarantees correctness. We use these parameters in the detailed description below. 
We now describe the cycles one by one. 

cycle 0-2: Same as cycles 0-2 in Section B.5.2. 

cycle 3: In the second pipeline stage, a carry-save representation of the the product N' Q <— NL\ * 
FLi is computed. In the first and third pipeline stages preparations are made for computing 
the product D[ <— D f Q • Fq in the next cycle. These preparations include the compression of 
D f Q to binary, the computation of the 3x multiple of D' Q , and the recoding of Fq. 

Note that the computation Fq «— 2 - {D' 0 [0 : 29] + 2~ 28 is done simply by complementing all 
the bits in the carry-save representation of D' Q [0 : 29]. The carry-save number representing 
Fq is then directly recoded without having to compress it. 

cycle 4: In the second pipeline stage, the product D[ <— D' 0 • Fq is computed. In the first and third 
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pipeline stages, preparations are made for computing the product N[ <— Nq • Fq in the next 
cycle. Namely, Nq is compressed and the 3x multiple of Nq is computed. 

cycle 5: In the second pipeline stage, the product N[ <— Nq • Fq is computed. In the third pipeline 
stage, D{ is compressed to reduce the error f x associated with the computation of F[. 

cycle 6: In cycle 6 preparations are made for computing N 2 <— N[ - F[ in the next cycle. These 
preparations include compressing N[> computing the x3 multiple of N[. As for F{, we 
rely here on the assumption that F{ e (1,1 + S x ] (see Section B.7.5). This implies that 
we only need to compute the lower bits of F[ (i.e., bits in positions greater than or equal 
to log 2 (l/<5i)). In Section B.6.2, we show that F{[1 : 26] is all zeros. Hence, we compute 
F{[27 : 56] simply by complementing the bits of D'[27 : 56]. We then proceed by recoding 
F{[27 : 56]. 

cycle 7: In the second pipeline stage, the product N% <— N[ • F[ + IN J rnd {mode) is computed. Our 
implementation is based on the fact that F[ is slightly larger than 1. Therefore we use the 
cheaper method of computing N[ + N[ • (F( - 1) instead of N x • F{. (See Section B.7.5.) 

In the first pipeline stage, preparations are made for the back-multiplication B • dewpoint — 
A. The definition of the dewpoint depends on N f 2 , which is not ready at this cycle. We there- 
fore use an approximation of the the dewpoint. We denote this approximation by dewpoint h . 
We derive dewpoint* 1 from the upper half of N[. (See also Section D.5.3 for more details.) 

In the third pipeline stage, we compute the two's complement of dewpoint h . We need this 
number for the next cycle for subtracting dewpoint* 1 . 

cycle 8: The first part of back-multiplication takes place in the second pipeline stage. We compute 
a H <— B • dewpoint h — A'. Note that, the final decision is based on positions [51 : 105] 
(see Section D.5.3). Hence a shortened representation of 4 — A suffices. Namely, instead of 
adding 4 - A\ we could simply add 2 50 - A 7 [51 : 52]. 
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In the third pipeline stage, we compute RC\ by compressing N% and truncating it at position 
52. 

In the first pipeline stage we compute the difference between the correct dewpoint and the 
estimate dewpoint* 1 . We denote this difference by dewpoint*. Note that a borrow-save 
representation of dewpoint* suffices since it only needs to be recoded in the first pipeline 
stage. See Section D.5.3 for more details on the computation of dewpoint*. We then recode 
dewpoint*, and the 3x multiple of B is prepared for the next cycle. 

cycle 9: The second part of the back-multiplication takes place in the second pipeline stage. 

cycle 10: In the third pipeline stage, (i) a compressed representation of RC 2 is computed by in- 
crementing RCu (ii) the 2-sign of the back-multiplication is determined from a L [51 : 105], 
and (ii) the IEEE rounded quotient is selected from RC X and RC 2 according to the rounding 
mode and the z-sign of a. 

B.6.2 Correctness 

In this section we prove the correctness of a design based on a rectangular 62 x 30 multiplier and 
an initial approximation with a relative error e 0 < 2~ 13,74 . To support double precision (pd = 53), 
dewpoint rounding requires that Assumption 19 holds. Namely, 0 < p^N^) < 2~ 53 . 

The analysis presented in this section is a concrete analysis for double precision. More gener- 
ally, this analysis applies to the following situation: (i) the accuracy of the initial approximation is 
roughly 2" p/2 ; (ii) the quotient is approximated by N' 2 (namely, the algorithm has 2 iterations of 
quadratical convergence); and (iii) the multiplier is a rectangular multiplier whose dimensions are 
not smaller than roughly p d x p<i/2- 

Claim 23 If the relative error in initial approximation satisfies |e 0 | < 2~ 13 74 , then correct double 
precision IEEE rounding is obtained with a rectangular 62 x 30-bit multiplier. 
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Proof: All we need to show is that 0 < piN^) < 2" 53 . By Theorem 15, 0 < p{N' 2 ) < tt 2 + 5 2 . 
We expand 7r 2 and 5 2 to obtain: 



To prove 7r 2 + S 2 < 2~ 53 , we prove upper bounds on the relative error terms (i.e., n i3 di) and 
substitute these upper bounds in Equation B.12. We bound these relative errors using the absolute 
error terms nepsi, deps^ fepst defined in Definition 6. 
We now bound each of the relative error terms: 

Bound on do: According to Table B.3, D' 0 is represented by a carry-save encoded digit string with 
digits in positions [0 : 60]. This digit string holds the rounded up value of D'_ Y • F'_ v This 
implies that the absolute error deps 0 associated with the computation of D' 0 is in the range 
[0, 2 • 2~ 60 ). Since the precise product D'_ x • FL X equals 1 — e 0 , this product is in the range 
1 ± 2- 13 - 74 . It follows that do < 2~ 59 /(l - 2" 13 74 ) < 2" 58 " 989 . 

Bound on n 0 : As in the case of D' 0 , Nq is also represented by a carry-save encoded digit string 
with digits in positions [0 : 60]. This digit string holds the rounded down value of N'_ x F!_ v 
This implies that the absolute error neps 0 associated with the computation of D' Q is in the 
range [0, 2 • 2~ 60 ). The precise product N'_i ■ FL X equals ~ • (1 - e 0 ), and is therefore at least 
(1 - |eo|). We conclude that n 0 < 2" 59 /(l - |e 0 |) < 2- 58 " 989 . 

Bound on / 0 : We first show that feps 0 < 2" 28 . The reason is that instead of computing 2 - 
jy o [0 : 60], we compute 2 - D' 0 [Q : 29] - 2" 28 . The absolute error is 2~ 28 - D' o [30 : 60]. 
Since the tail 0 < D' o [30 : 60] < 2 • 2~ 29 , we conclude that indeed 0 < feps 0 < 2" 28 
(recall that D'[0 : 60] is represented by a carry-save encoded digit string). Next, we expand 



1 - (1 - n 2 ) • 



1 — n 0 1 — rii 



7T2 



1 + d 0 * 1 + d x 




(B.12) 
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2 - D' 0 = 2 - (1 + d 0 ) • (1 - e 0 ). We conclude that 

/o ~ 2 -(l + do)-(l + |e 0 |) 

< 2-27.999894 



Bound on d\i As in the bound on d 0 , depsi < 2 59 . We bound from below the exact product 
as follows: D' 0 ■ F£ = (1 - / 0 ) • D' 0 • (2 - Now Z% = (1 + d 0 )(l - ^o), hence 

#o e [1 - |e 0 |, 1 + S 0 ). It follows that D' Q • > (1 - / 0 ) - (1 — 5g) > 1 — <*i. We conclude 
that 

2-59 
< 2-58.999999 



Bound on rt\\ Again, we use nepsi < 2~ 59 . We bound from below and F$ as follows: Nq 
(1 - n 0 ) • (1 - co) • Hence, iVJ > (1 - n 0 ) ■ (1 - |e 0 |). Also, Fq — (1 — f 0 ) • (2 
£>o) > (1 - /o) • (1 - <*o). We conclude that 



2~59 

Tli < 



(l-no).(l-|co|)-(l-/o)-(l-*) 

< 2~ 58.999789 



Bound on fx: Assume that D[ e [1 — 5i, 1). Note that according to Table B.3, when computing 
F[ (during cycle 6) we have a binary non-redundant representation of D[ (from cycle 5). 
Since 

Si = &l + fo 

< 2~26.716654 
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It follows that the binary representation of D[ satisfies D[ [0 : 26] = 1 - 2" 26 . We conclude 
that 

2 - D'[0 : 60] = (l + (1 - 2~ 26 ) + 2~ 26 ) - (D[[0 : 26] + D[[27 : 60]) 
= 1 + 2" 26 -D[[27 : 60], 

Instead of computing 2 - £>J [0 : 60] , we compute (in cycle 6) 1 + 2~ 26 - (D[ [27 : 56] + 2" 56 ) 
(note that the 1 does not appear explicitly since it does not effect the actual computation). It 
follows that 

fepsx = (1 + 2" 26 - D[[27 : 60]) - (l + 2~ 26 - (D[[27 : 56] + 2" 56 )) 
= 2" 56 - &[57 : 60] 
< 2 -56 . 

Since we assumed that D[ < 1, it follows that 2 - D[ > 1. Hence /i < fepsi < 2 -56 . 
If D\ > 1, then, by Equation B.3 in the proof of Theorem 15, we know that 

< 2~56.999894 

However, by Assumption 14 (precise multiplication by 1), it follows that N[ = Hence, 
if D[ > 1, then we obtain the desired relative error. 

Bound on ra 2 : As before, neps 2 < 2~ 59 . 

The precise product N[ • F[ is bounded as follows. Saturated rounding implies that F[ > 1. 
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Theorem 15 implies that N[ > ^ • (1 - 7Ti - <5i). We conclude that 

2~59 
(1 - 7Ti - Si) 
< 2-57.9999999 

We now combine these bounds to obtain: 

tt 2 + 8 2 < 2- 5305992 < 2~ 53 . 

□ 

We point out that Claim 23 is tight in the sense that if we decrease any of the three parameters: 
initial approximation error, long operand length, or short operand length, then we can no longer 
prove that p(N' 2 ) < 2 -53 . 

B.7 Hardware Optimizations 

B.7.1 Redundant Feedback & Partial Compression 

Addition trees in multipliers are not amenable to pipelining, so short clock cycles are not achievable 
with reasonable cost if the addition tree has too many rows. Booth radix 8 recoding reduces the 
number of rows in an addition tree from n to • 

Booth radix-8 multipliers are usually implemented using a 3-stage pipeline (i.e. (i) precom- 
pute the 3x multiple of the first operand of the multiplier and recode the second operand, (ii) 
addition tree, and (iii) final carry-propagate addition). Goldschmidt's algorithm performs only 2 
multiplications per iteration. Hence running Goldschmidt's algorithm on a 3-stage pipeline cre- 
ates un-utilized cycles (i.e. "bubbles") in the pipeline. The Athlon™ processor (in hardware) [28] 
and Itanium™ processor (in software) [5] allow other multiplications to be executed during such 
bubbles. 
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Following Seidel [35], we manage to get rid of all of these bubbles but one by allowing redun- 
dant operands. Conceptually, we skip the third pipeline stage and reduce the pipeline to a two-stage 
pipeline for all but the last iteration of the algorithm. This is obtained by feeding back results rep- 
resented in carry-save format. This design is not symmetric in the sense that the first operand and 
second operand of the multiplier are processed differently during the first pipeline stage. During 
the first pipeline stage, operands represented as carry-save number are processed as follows: 

1. The first operand is compressed and its 3x multiple is computed. This requires two adders. 
To save hardware, we suggest "borrowing" the adder from the third pipeline stage for the 
purpose of compressing the first operand. This explains why we listed the compression of 
the first operand as an operation that takes place in the third pipeline stage in Tables B.2 
and B.3. 

2. The second operand can be partially compressed from carry-save representation before being 
fed to the Booth recoder. 

B.7.2 Normalization by checking (A > B) 

Suppose that the operands satisfy \A\, \B\ e [1,2). We compare the operands A and B. If \A\ < 
\B\ y then we set \A'\ *- 2 • \A\; otherwise \A'\ <- This ensures that \A'\/\B\ e [1,2) is the 
value of the normalized quotient. This normalization also slightly improves the bound on the error. 

B.7.3 Directed rounding of carry-save numbers 

Our algorithm has to deal with rounding in two situations. The first situation is when directed 
rounding is used in intermediate computations. In this situation, the intermediate product is repre- 
sented by a carry-save encoded digit string. We wish to apply directed rounding without incurring 
a significant delay. Directed rounding of carry-save number is described in this section. The sec- 
ond situation is when we wish to determine the correct IEEE rounded quotient; we deal with IEEE 
rounding in Section B.4. 
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We propose an implementation of Algorithm 2 in which intermediate results are represented 
by carry-save encoded digit strings. Since Algorithm 2 uses directed roundings for all intermediate 
computations, we need to explain how directed rounding is applied to carry-save numbers. 

Notation. Consider three binary vectors x, y, z. We denote 3:2-addition by FA(x, y, z). Namely, 
FA(x, y y z) means that the three binary vectors x, y, z are fed to a line of full-adders. The output 
of FA(x, y, z) is two binary vectors s and c (i.e., a carry-save number) that satisfy \x\ 4- \y\ + 

N = N + |c|. 

By truncating a carry- save encoded digit string a[i : j] after position q we simply mean chop- 
ping off the tail and leaving only a[i : q]. We denote truncation of a after position q by l<r\ g . 

We often regard a carry-save encoded digit string a also as two binary vectors a' and cr". 
Hence if a is a carry-save encoded digit string and a: is a binary vector, then FA(a, x) simply 
means FA{a\ a", x). 

The following lemma deals with rounding up of a carry-save number. It shows that 3:2-addition 
followed by truncation after position q can be used to implement approximate rounding up. The 
3:2-addition adds a constant, called the injection, that depends only on the rounding position and 
the length of the number. 

Lemma 24 Let a[0 : t] denote a carry-save encoded digit string. Let INJCS £) E [2~ 9 + 
2-(«+ 1 ) ) 2 • (2~ q — 2~*)] and assume that W}CS RU (q, t) is represented by a binary string I[q : t). 
Then, 

\a[0:t}\ < \[FA(a[0:t]J[q:t])} q \ < \a[0 : q}\ + 2~* +1 . 
Proof: We first show that 



\a[0:t}\ < \[FA(a[0:t],I[q:t])\ q \. 



(B.13) 
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From the range of iNJCSncfa, t), it follows that I[q] = 1 and I[q + 1] = 1. Hence, 

\[FA(*[0 : t], I[q : t])) q \ = \ a [0 : q]\ + 2"* + C ■ 2"*, (B.14) 

where C denotes the carry bit generated by adding /[? + 1] + <j[q + 1]. We consider two cases: (i) 
C = 0: in this case a[q + 1] = 0, and hence, \a[q + 1 : *]| < 2^. It follows that Equation D.9 
holds in this case, (ii) C = 1: in this case a[q + 1 : t] < 2~^~ x \ and hence Equation D.9 also 
holds. 

We now show that 

\IFA(*[0 : t)J[q : t])\ q \ < \a[0 : q]\ + 2~*+\ (B.15) 

This follows directly from Equation D.10 since \a[0 : q]\ + + C • 2" 9 < |<j[0 : $]| + 2~ q+1 . □ 

We present in Lemma 43 the widest possible interval for INJCS ^(q, t). It is probably most bene- 
ficial to use a value with the smallest number of ones, namely, iNJCS/^(g 5 1) = 2~ q + 2~ (9+1) . 

B.7.4 Implementation of Saturated Rounding 

The definition of saturated rounding from Assumption 2 requires that F[ > 1 for i > 0. Suppose 
that D$>1, then 2 - D\ < 1, and hence F[ should equal 1. Moreover, Assumption 14 requires 
that every multiplication by 1 is precise. 

We propose the following implementation. Multiplication takes place in the second pipeline 
stage. During a multiplication by F[ that is part of iteration i > 0, the algorithm checks if D\ > 1. 
If A' ^ 1. then we simply de-activate the clock enable signal of the register that is supposed 
to store the result of the multiplication. The effect of ignoring the multiplication means that the 
multiplicand is not changed, i.e., precise multiplication by 1 is obtained. 

We test whether D[ > 1 as follows. Recall that is represented in carry-save format. We 
subtract 1 from the D\ using a signed 3:2-adder (i.e., PPM-add [26]) to obtain a borrow-save 
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number. The borrow-save number is fed to a subtracter, and the decision is made according to the 
value of the sign bit of the difference. 

B.7.5 Half-size multiplication in the last iteration 

The precision in the last iteration must be at least p d to obtain a sufficiently close approximation 
of the quotient. This presumably implies that the length of the second operand (i.e. the multi- 
plier) should be roughly p d . In Section B.6 we present an implementation that only uses a single 
pass through a half-sized multiplier. This implementation is based on the error analysis and, in 
particular, on Corollary 17. 

For the sake of concreteness we consider double precision (p d = 53). The last iteration requires 
computing the product N 2 = N[ • F{. The length of F{ should be bigger than p d = 53 by a few 
bits; otherwise the relative error f{ would be too large. How can we avoid having to multiply N[ 
by the full length of F{7 

Rewrite N 2 as follows: 

N' 2 = N[ • (1 + (F{ - 1)) (B.16) 
= N[ + N['(F{-1). (BM) 

By Corollary 17, (F[ — 1) € [0, Si]. Hence, we may ignore bit positions [1 : j] where j = |k>g 2 ^-J • 
We conclude that N 2 = N[ + N[ • F[\j : t], where t = flog 2 Our error analysis shows that 
t — j is roughly p d i% as required. 

The addition of N[ to the product N[ • F{[j : t] can be done by using one of the extra rows 
in the adder tree. The second extra row is used for feeding the rounding injection lNJ rnd (mode). 
Note that RZ is used for the computation of N 2 , hence no carry-save injection is required. 
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B.7.6 Optimized implementation of Dewpoint Rounding and Back multipli- 
cation 

In Section B.4.2 we showed how the IEEE rounded quotient is chosen to be either RC X or RC 2 
based on the z-sign (positive, negative, zero) of dewpoint • \B\ — \A'\. 

For the consideration of negative numbers and signs, we need to consider two's complement 
representations. In particular we consider a representation called two's complement carry-save 
representation [7]. In this representation, a number is represented by two binary vectors, each of 
which is a two's complement number, namely, the most significant position has a negative weight. 
Our goal in the back multiplication is to compute the sign of a = \B\ ■ dewpoint — \A'\, we also 
want to know if a is precisely zero to determine the condition for rounding mode RI. The z-sign of 
a number x is either positive, negative, or zero depending on the result of comparing x with zero. 
Our goal is to compute z-sign(a). 

We now discuss two optimization techniques. The first technique is used to simplify the 
computation of z-sign(a). The second technique is used in double precision to split the back- 
multiplication into two passes through the half-sized multiplier, where in the first pass an estimate 
of the dewpoint is used. The estimated dewpoint is based on N[ rather than N' 2 . This means that 
we begin the first pass of the back-multiplication before N£ is ready. 

Computing z-sign(a). We first prove the following claim. 
Claim 25 

-2 - 2~ p < \B\ • dewpoint - \A'\ < 2 • 2" p . 
Proof: By Assumption 19, 0 < p(N' 2 ) < 2~ p . Recall that 



dewpoint = [iV^ + lNJ rmi (mode)Jp_i + C^imode). 
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It follows that 

dewpoint < N 2 + INJ rnd (mode) + C^mode) 

< N' 2 + 2" p+1 
A 1 

<4r + 2- p+1 . 

The second line follows from the definition of the rounding injection inj^ (mode) and the dew- 
point displacement constant C^imode), The third line follows from 0 < p(N^). 

Since p{N f 2 ) < 2"* and A'/B < 2, the other part of the claim follows if we show that 
dewpoint > N 2 . We consider each of the three rounding modes. 

RZ: Let x and x' denote two successive represen table significands that sandwich N 2 , namely, 
N' 2 e [x,x') where x' = x + 2~p +1 . In RZ, \W rnd (RZ) = 0 and C^ W {RZ) = 2"P +1 . 
Hence: [N 2 + iNJ rnd (mode)J p _i = x and dewpoint = x*. 

RI: Assume here that N£ e (x, x f ] (i.e., the equality can hold in x' instead of x). In RI, 
lNJ rnd (i?/) = 2"P +1 - 2~ l and C^RNU) = 0. Hence [N£ + lNJ rnd (mode) J p _! = x' and 
dewpoint > x r . 

RNU: Assume here that e [x + 2"", x' + 2"P) (i.e., we sandwich ^ between two odd 
multiples of 2-P In RNU, iNJ rnd (,RiVi7) = 2~p and C^RNU) = 2"*\ Hence [JV* + 
lNJrnd(7node)J p _ 1 = x' and dewpoint = x' + 2" p . 

In all three rounding modes, dewpoint > N 2 , and hence the claim follows. □ 

Let a[— 2 : £] denote a two's complement non-redundant representation of a, namely: 

£ 

a = -a_ 2 • 2 2 + ^ a, • 2~\ 
t=-i 
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By Claim 44 it follows that -a_ 2 -2 2 +£?J?i a i* 2~* e {-2 p - 2 ,0}. Hence, the bits in positions 
[—2 : p — 2] are either all zeros or all ones. We conclude that 

a = 0 iff OR(a p _2,...,at) = 0, (B.18) 
a<0 iff a p _ 2 = 1. (B.19) 

Note that the multiplier computes a two's complement carry save representation of a. For the 
purpose of computing ^-sign(a), we conclude that the carry-save digits of a in positions [p — 2 : t] 
suffice. This last remark also implies that we need not compute 4 — A! since it suffices to use 
2-(p-3) _ A >fc _ 2 : p - l]. 

Back-multiplication with an estimated dewpoint. In Section B.6, back-multiplication is im- 
plemented by two passes through the half-sized multiplier. An estimate of the dewpoint is used in 
the first multiplication. We denote this estimate by dewpoint* 1 . The estimate dewpoint h is com- 
puted by truncating N[ to fit the length of the second operand of the multiplier. In the first pass we 
compute 

a H <- B • dewpoint* 1 + (4 - A'). (B.20) 

In the second pass we correct the computation. The difference between the dewpoint and the 
estimated dewpoint is denoted by dewpoint 1 . Formally, 

dewpoint 1 = dewpoint — dewpoint h 

= RCi + Cde W (mode) — dewpoint h . 

In the second pass we compute 



a L <— B • dewpoint* + a H . 



(B.21) 



Note that Equations B.20 and B.21 imply that a* = a. Based on Equations B. 18 and B. 19, we are 
only interested in a L \p - 2 : t]. 

We now compute which digit positions are required for dewpoint'. To be concrete we focus 
on double precision as analyzed in Section B.6. The estimate dewpoint* is simply the truncation 
of the non-redundant representation of N{. Since the multiplier's second operand length is 30, we 
set 

dewpaint h <— N[[0 : 29]. 



Claim 26 Under the premises of Claim 23, 0 < dewpoint 1 < 2~ 25A47 . 

Proof: The rounding candidate BC X equals the truncation of A*[ 0 : 60] + lw rnd (mode) at 
position 52 (note that the carry bit for position [52] resulting from compressing the carry-save 
representation of JVJ is computed). Note that A* = (1 - n 2 ) • F{ • N[ . Since F[ < 1+S U it follows 
that^ < (l + s^-Nl Hence, 

dewpoint? < (1 + S x ) ■ N[[0 : 60]| + !NJ rnd (mod e ) + C^mode) - N'[0 : 29] 
< (1 + St) ■ N[[30 : 60] + J, • N[[0 : 29] + 2~ 52 . 

The last equation follows from iW rnd (mode) + C^mode) < 2~^\ We bound N[[0 : 29] < 2, 
NtfO : 60] < 2-» and 6> < 2~ 2 ^ (from me proof of C]aim 23) tQ ^ ^ ^ J 
dewpoint'. The lower bound follows simply because #d > dewpoint h . D 
Claim 26 implies that if dewpoint 1 is represented as a carry-save number, then it suffices to 
represent dewpoint 1 using digits in positions [26 : 53]. This means that dewpoint 1 easily fits as the 
second operand of the multiplier (if shifted by 26 positions). Finally, since we are only interested 
in positions to the left of p d - 2, we conclude that Equations B.20 and B.21 can be implemented 
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as follows: 

a H [0 : 81] <- B[0 : 52] • dewpoint h [0 : 29] + (2 50 - A'[51 : 52]) 
a L [26 : 105] <- B[0 : 52] ♦ dewpoint £ [26 : 53] + a" [51 : 81]. 

Note that this means that in the second pass we need to also shift a H [51 : 81] by 26 positions 
before we inject it to the adder tree. 

Computation of dewpoint*[26 : 53] The implementation of the dewpoint computation is 
depicted in figure B.3. We explain the details of the implementation below. The operand 
dewpoint*[26 : 53] is required in cycle 8 for the double precision implementation (see table B.3). 
A Booth recoded representation of the value dewpoint* [26 : 53] needs to be computed from the 
carry-save representation of N% [0 : 60]. The computation is based on the equation: 

dewpoint e [0 : 53] = A^fO : 52] + carry 52 ( N%[53 : 60]) + Cde W (mode) - dewpoint h [0 : 29]. 

By Claim 26, 0 < dewpoint* < 2" 25 . Hence, it suffices to compute bit positions [26 : 53] of the 
binary representation of dewpoint*. However, computing the binary representation is slow and is 
not needed for obtaining a Booth radix-8 recoding of dewpoint*. 

We precompute (part of) 2 — dewpoint h [0 : 29] in cycle 7. Consider a carry-save number 
RCS[-1 : 53] defined by 

RCS[-1 : 53] = N%[0 : 52] + carry 52 (N%[53 : 60]) + C^mode) + (2 - dewpoint h [0 : 29]). 

Obviously, the value represented by RCS[— 1 : 53] equals 2 + dewpoint*. 

Our goal is to compute a Booth radix-8 recoding of dewpoint*. We denote a Booth radix-8 
recoding of dewpoint* by BD^l— 1 : 53], where BDd^lZi + 2] e [—4,4] is a Booth radix-8 
digit, for every i = —1,0,. ..,17. First, we observe that (i) BD^l— 1 : 23] = 0, namely the 
Booth digits in positions [—1 : 23] are all zeros, and (ii) the most significand non-zero Booth digit 
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in BDdeul— 1 : 53] is positive. This follows from the fact that dewpoint* 6 [0, 2~ 25 ) and since 
the fraction range of a Booth radix 8 number is (—4/7, 4/7). Hence, we are left with the task of 
computing BD^ [26 : 53]. 

We apply the following method of [35] for obtaining a minimally redundant format in which 
each block of 3 positions represents a value in the range [—4,4]. First we feed the carry-save 
number RCS to a line of 3-bit adders. The output of the line of 3-bit adders for positions [3z : 3i + 
2] consists of 4 bits: one bit for position [3i] has weight 2" 3i = 4 • 2~< 3i+2 \ one bit for position 
[3i+ 1] has weight 2~ 3i ~ l = 2-2"( 3i+2 > , and two bits for position [3z+2] each have weight 2~< 3i+2 ). 
It follows that this 4-bit string represents a digit in the range [0, 8] whose weight is 2~ (3i+2) . To 
allow for interpretation of each string as a digit in the range [-4, 4], we add to RCS a constant. 
Loosely speaking, the constant increases each block [3i : 3i + 2] by 4 • 2~< 3i ~ 2 >. This implies that 
subtracting 4 from each digit is in effect a subtraction of the added constant. By subtracting 4 from 
each digit, we translate digits in [0, 8] to Booth radix-8 digits in [-4, 4]. We denote the Booth digit 
correspoding to block [3i : 3i + 2] by BD[3i + 2]. 

Based on this discussion our circuit does not compute RCS, Instead, the algorithm computes 
(part of) a carry-save number R[— 1 : 53] that is defined by 

R[-l : 53] = RCS[-1 : 53] + 2 3 + 2° + • • ■ + 2" 24 + 2~ 27 + • • • + 2" 51 . 

The carry-save number i?[24 : 53] is fed to a line of 3-bit adders that outputs (part of) the sum- 
string Sz[— 1 : 53] and (part of) the carry-string Cs[— 2 : 52]. Note that the carry-string C3 has bits 
only in position 3i + 2, for each block [3i : 3i + 2]. 
We define the Booth digit BD[3i + 2] by: 

BD[3i + 2] = S 3 [3i : 3i + 2] + C 3 [3i + 2] - 4. 

Note that the above discussion in which a constant is added and then subtracted from RCS implies 



96 



that RCS[-1 : 53] = BD[-1 : 53]. 

Finally, we define BD'[26] as the Booth digit in position [26] of R! = R - 2~ 24 . Of course, 
i?D'[26] is not uniquely defined, so we describe how we compute it. The Booth digit BD* [26] is 
computed by not adding the constant 2" 24 in the left most 3-bit adder (see Fig. B.3). We then use 
a 3-bit adder to compute the sum bits S3 [24 : 26], and define 

BD'[26] = S' s [24 : 26] + C 3 [26] - 4. 

Claim 27 Let £[24 : 29] denote the binary representation ofN%[25 : 30] + 2 -24 - dewpcrint h [25 : 
29]. Then, 



dewpoint e = BD deiu [26 : 53] = < 



BD[26 : 53] ifB[24 : 25] = 00 

BD'[26) o BD[29 : 53] otherwise. 



Proof: To shorten notation we define a, /3, and 7 as follows: 

a = N%[-1 : 24] + 2 - 2~ 24 - dewp<rint h [Q : 24] 

P = N$[2S : 30] + 2~ 24 - dewpoint h [25 : 29] 

7 4 A£[31 : 52] + carry 52 (N%[53 : 60]) + C^mofe). 

Note that <* + /3 + 7 = 2 + dewpoznt* and that /? + 7 = £L>[26 : 53]. 

1. B[24 : 25] = 00. In this case p < 2~ 25 - 2" 29 . Since 7 < 2" 29 + 2" 53 , it follows that 
P + 7 < 2 -24 . Now a is a multiple of 2 -24 , and 2 4- 2 -25 >a + /3 + 7>2. Hence a = 2. 
It follows that dewpoint e = /? + 7 = Z?Z)[26 : 53], as required. 

2. £[24 : 25] = 01 or £[24 : 25] = 10. In this case 2~ 25 < p < 2~ 24 + 2" 25 . It follows 
that 2~ 25 < P + 7 < 2 • 2 -24 . If a > 2, then a + /? + 7>2 + 2~ 25 , a contradiction. If 
a < 2 - 2 • 2~ 24 , then a + P+ 7 < 2, a contradiction. Hence, a = 2 - 2~ 24 . We conclude that 
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dewpoint* = a+/3+ 7 -2 = /?+7-2~ 24 . Note that BD'[26]oBD[29 : 53] = a-f/3-2" 24 , 
as required in this case. (There is a subtle issue which we should address here as well as in 
the previous case. Namely, show that C3[23] = 0. This follows from dewpoint e G [0, 2" 25 ). 
If C£[23] = 1, then dewpoint* > 2" 23 • (1 - 4/7) > 2~ 25 , a contradiction.) 

3. £[24 : 25] = 11. In this case 2~ 24 + 2" 25 < (3 < 2" 23 - 2" 28 . Therefore, 2" 24 + 2~ 25 < (3 + 
7 < 2~ 23 . If a > 2 - 2 -24 then a + /? + 7>2-l- 2" 25 , a contradiction. If a < 2 - 2" 23 then 
a + /? + 7<2, a contradiction. Hence, this case cannot occur. 

□ 

Finally, note the same hardware can also be used for single precision implementation by unify- 
ing and aligning the rounding position for both cases with a 3-bit right shift in the single precision 
case. The implementation depicted in figure B.3 and figure B.2 presents such an integration for 
single and double precision. 

B.8 Pipelining options 

In Section B.6 we presented a micro-architecture that implements our algorithm. The core of this 
three stage micro-architecture is a half-sized multiplier. 

In Tables B.2-B.3 we showed that single precision division can be computed in 9 clock cycles 
and double precision division can be computed in 11 clock cycles (this includes one cycle for 
reciprocal approximation). It is easy to see that one can schedule a new independent division every 
6 cycles in single precision and every 8 cycles in double precision. In Table B.l, we refer to this 
micro-architecture as "proposed design (1/2)". 

We now discuss how to reduce the restart-times between independent divisions; this requires 
additional hardware to handle the additional amount of work. 

In part (A) of Figure B.4 we present a time-space diagram for double precision division (i.e., 
p d = 53) based on the schedule defined in Table B.3. The x-axis refers to clock cycles and the 
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y-axis refers to the four pipeline stages. The squares of the diagram are labeled by the smallest 
hardware that can execute the corresponding task: A l/4 denotes reciprocal approximation that 
returns roughly p d /4 bits of precision (as well as a comparator), Mj /2 denotes the ith pipeline stage 
of a half-sized multiplier, and M[ /A denotes the ith pipeline stage of a quarter-sized multiplier. 

Reduced restart-times can be obtained by using more than one multiplier. We propose three 
designs nicknamed Proposed Designs (1), (2), and (3). Note that a rough cost estimation shows 
that these designs require about 1 x , 2 x , and 3 x the hardware cost of a conventional multiplicative 
division implementation (whose core is full-sized multiplier). 

Proposed Design (1) Proposed design (1) is built of two copies of half-sized multipliers. Part (B) 
of Figure B.4 depicts the assignment of tasks to the two multipliers. We refer to the first multiplier 
by Ml and to the second multiplier by M2. The restart-time of this design is only 4 cycles since 
Ml and M2 are engaged in the computation for four cycles and three cycles, respectively. 

Proposed Design (2) Proposed design (2) is built of three copies of half-sized multipliers and 
one copy of a quarter-sized multiplier. Part (C) of Figure B.4 depicts the assignment of tasks to 
the four multipliers. Here we refer to the quarter-sized multiplier by Ml and to the three half- 
sized multipliers by M2, M3, and M4. The restart-time of this design is only 2 cycles since each 
multiplier is engaged in the computation for at most two cycles. 

Proposed Design (3) Proposed design (3) is built of five copies of a half-sized multiplier and two 
copies of a quarter-sized multiplier. Part (D) of Figure B.4 depicts the assignment of tasks to the 
four multipliers. Here we refer to the quarter-sized multipliers by M1-M2 and to the five half-sized 
multipliers by M3-M7. This design is fully pipelined, and allows for a new independent division 
every cycle. 

The latencies and restart times of the four proposed designs are listed in Table B.l. These 
designs provide trade-offs between hardware cost and division throughput. The cost ranges from 
roughly one half up to three times the cost of a regular multiplicative division implementation. The 
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throughput for double precision ranges from one division every 8 cycles up to the fully pipelined 
implementation that allows to start a new division every cycle. 

Similar reductions in restart-times can be applied to the micro-architecture presented in Sec- 
tion D.6. We focus on the situation in which a roughly p d x multiplier is used both for 
double-precision (i.e., p d = 53) and single precision (i.e., p 3 = 24). When performing single 
precision division, the multiplier is a full-sized multiplier. (If fact, the multiplier is two full-sized 
multipliers.) By applying the same method to the schedule given in Table B.2, the restart-times 
reported in Table B. 1 are obtained. 
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| cycle | RTL 


1 comment 


0 


R2<r-approx(l/B) 






Rl *— B 




1 


R7 <- recode(R2) 






/to < — o - rtl 


= 3B 






— b 




#4, m <- 2,4 if A < otherwise A 


= A' 


2 


#6 <- 3 • m 


= 3 A' 




i?5 <- Rl 


= A' 




R\ <r- B 






R8<r-R5 R7 + iNJ Ht/ (60, 89) 


= D'o 


3 


R3 <— compress(RS) 


= D' 0 




R8+-R5- R7 


= K 


4 


Rb <— compress(R8) 


= Atf 




#6 <- 3 ■ #8 


= 3*5 




R7 recode( two's complement(RZ)) 




5 


R8 *— R5 ■ R7 -f iNj rrid (mode) 


= N{ 


6 


727 <— recode(shifi s (dewpoint(R8))) 


dewpoint 




R5<- Rl 


= B 




R6<^3Rl 


= 3B 




RIO <— compress(R8) 


= RCi 


7 


R8<r- R5 R7+ {two's complement(R4)[22 : 23]) 


= a 


8 


RC 2 = RIO + 2~ 2a 


compressed RCi is ready 




*est(i*8[22 : 47]) 


test(a) 



Table B.6: Register transfer language describing execution of a single precision division in micro- 
architecture described in Figure B.2 
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| cycle | RTL | comment ~| 



o 


R2 <— armroxfl IB\ 






Rl B 




1 

± 


ill ~ / tZLAJUt\ZyJ.V£i J 






R6 <— 3 • Rl 

J. *,\J ^ O ill 


= 3B 




il5 <— ill 


= B 




il4, ill <— 2/1 if A < J9, otherwise A 


— A 1 


2 


il6 «— 3 • ill 


= 3- A' 




il5 <— ill 


= A' 




ill <— £ 






il8 «- il5 • il7 + INJhc/(60, 89) 


= D' 0 


3 


R7 * — recode(ones-complement(R8)) 


= FA 

u 




R6 <— 3 • il8 


= 3 D' 0 




il5 < — compress^ il8) 


= DL 
u 




il8 <— il5 • il7 


= M 

u 


4 


il5 « — CO?7ip7*6Ss(il8) 


= N{ S 




il6 <— 3 • il8 


= 3N{> 




il8 <— il5 • il7 + lNJflt/(60, 89) 


= D[ 


5 


il8 4— il5 • il7 


= N[ 




il3 < — compress ( R8) 


= D[ 




Rll <- (il8 < 1) 


for saturated rounding in cycle 7 


6 


il7 <— recode(shift(ones-complement{R3))) 


F{[27:56\ 




il6 <- 3 • il8 


= 3N{ 




il3, il5, illO <- compress(il8) 


= N{ I 


7 


il7 recode(il3) 


= dewpoint* 1 




il6 3 ill 


= 3- B 




il5 <- ill 


= B 




il3 2 - illO 


prepare (2 — dewpoint*) 




if illl then il8 <- il5 • il7 + shift(Rb) + lNJ rn d(mo<fe) 


= N£ 


8 


illO <- compress(il8) 


ilCi compressiNg) 




il7 <— dewpoint-recode(R8, il3, C^ew) 


compute & recode dewpoint c 




il8 4- il5 • il7 + (2 - il4) 


= a",{2- R4) = 2- A' 


9 


il8 <- il5 • il7 + il8[51 : 81] 


= a^ 


10 


RC 2 = RIO + 2~ M 


compressed ilC2 is ready 




test(il8[51 : 105J) 


testis ) 



Table B.7: Register transfer language describing execution of a double precision division in micro- 
architecture described in Figure B.2 
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stage 3 



rmode 




rounded quotient 



Figure B.2: Microarchitecture for dual mode single precision and double precision Division Im 
plementation. 
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Proposed design (3) 
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Figure B.4: (A) A time-space diagram of double precision division using the micro-architecture 
depicted in Fig. B.2. A x / 4 denotes reciprocal approximation with precision roughly p/4, M{ /2 
denotes the zth pipeline stage of a half-sized multiplier, and MJ /4 denotes the ith pipeline stage 
of a quarter-sized multiplier. Proposed design (1): Ml and M2 denote two half-sized multipli- 
ers. Proposed design (2): Ml denotes a quarter-sized multiplier, and M2-M4 denote half-sized 
multipliers. Proposed design (3): M1-M2 denote quarter-sized multipliers, and M3-M7 denote 
half-sized multipliers. 
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Appendix C 

A Parametric Error Analysis of 
Goldschmidt's Division Algorithm 

(Joint work with Warren E. Ferguson. An extended abstract appeared in the proceedings of the 
16th IEEE Symposium on Computer Arithmetic. Full version submitted to Journal of Computer 
and System Sciences.) 

abstract 

Back in the 60's Goldschmidt presented a variation of Newton-Raphson iterations for division that 
is well suited for pipelining. The problem in using Goldschmidt's division algorithm is to present 
an error analysis that enables one to save hardware by using just the right amount of precision for 
intermediate calculations while still providing correct rounding. Previous implementations relied 
on combining formal proof methods (that span thousands of lines) with millions of test vectors. 
These techniques yield correct designs but the analysis is hard to follow and is not quite tight. 

We present a simple parametric error analysis of Goldschmidt's division algorithm. This anal- 
ysis sheds more light on the effect of the different parameters on the error. In addition, we derive 
closed error formulae that allow to determine optimal parameter choices in four practical settings. 
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We apply our analysis to show that a few bits of precision can be saved in the floating-point 
division (FP-DIV) micro-architecture of the AMD-K7™ microprocessor. These reductions in pre- 
cision apply to the initial approximation and to the lengths of the multiplicands in the multiplier. 
When translated to cost, the reductions reflect a savings of 10.6% in the overall cost of the FP-DIV 
micro-architecture. 

C.l Introduction and Summary 

Asymptotically optimal division algorithms are based on multiplicative division methods [23, 30, 
37]. Current commercial processor designs employ a parallel multiplier for performing division 
and square-root operations for floating-point [1,5, 22, 28]. The parallel multiplier is used for 
additional operations, such as: multiplication and fused multiply-add. Since meeting the precision 
requirements of division operations requires more precision than other operations, the dimensions 
of the parallel multiplier are often determined by precision requirements of division operations. It 
follows that tighter analysis of the required multiplier dimensions for division operations can lead 
to improvements in cost, delay, and even power consumption. 

The main two methods used for multiplicative division are a variation of Newton's method [14, 
15] and a method introduced by Goldschmidt [16] that is based on an approximation of a series 
expansion. Division based on Newton's method has a quadratic convergence rate (i.e., the number 
of accurate bits doubles in each iteration) and is self-correcting (i.e., inaccuracies of intermediate 
computations do not accumulate). A rigorous error analysis of Newton's method appears in [3, 20, 
25] and for various exceptional cases in [5]. The analysis in [3, 20] considers the smallest precision 
required per iteration. Our error analysis follows this spirit by defining separate error parameters 
for every intermediate computation. In addition, the analysis in [3, 20] relies on directed roundings, 
a method that we use as well. 

Each iteration of Newton's method involves two dependent multiplications; namely, the prod- 
uct of the first multiplication is one of the operands of the second multiplication. The implication of 
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having to compute two dependent multiplications per iteration is that these multiplications cannot 
be parallelized or pipelined. 

Goldschmidt's division algorithm also requires two multiplication per iteration and the conver- 
gence rate is the same as Newton's method. However, the most important feature of Goldschmidt's 
algorithm is that the two multiplications per iteration are independent and can be pipelined or 
computed in parallel. On the other hand, Goldschmidt's algorithm is not self-correcting; namely, 
inaccuracies of intermediate computations accumulate and cause the computed result to drift away 
from the accurate quotient. Goldschmidt's division algorithm was used in the IBM System/360 
model 91 [2] and even more recently in the IBM S/390 [33] and in the AMD-K7™ microproces- 
sor [28]. However, lack of a general and simple error analysis of Goldschmidt's division algorithm 
has averted most designers from considering implementing Goldschmidt's algorithm. Thus most 
implementations of multiplicative division methods have been based on Newton's method in spite 
of the longer latency due to dependent multiplications in each iteration [5, 22] (see also [38] for 
more references). 

Goldschmidt's method is not self-correcting as explained in [21] (there is a wrong comment 
on this in [39]). This makes it particularly important and difficult to keep track of accumulated 
and propagated error terms during intermediate computations. We were not able to locate a gen- 
eral analysis of error bounds of Goldschmidt's algorithm in the literature. Goldschmidt's error 
analysis in [16] is with respect to a design that uses a serial radix-4 Booth multiplier with 61-bits. 
Goldschmidt's design computes the quotient of two binary numbers in the range [1/2, 1), and his 
analysis shows that the absolute error is in the range [— 2 56 , 0]. Krishnamurthy [21] analyzes the 
error only for the case that only one multiplicand is allowed to be imprecise in intermediate compu- 
tations (the second multiplicand must be precise); such an analysis is only useful for determining 
lower bounds for delay. Recent implementations of Goldschmidt's division algorithm still rely on 
an error analysis that over-estimates the accumulated error [28]. Such over-estimates lead to cor- 
rect designs but waste hardware and cause unnecessary delay (since the multiplier and the initial 
lookup table are too large). These over-estimations were based on informal arguments that were 
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confirmed by a mechanically verified proof that spans over 250 definitions and 3000 lemmas [31]. 

Agarwal et al [1] presented a multiplicative division algorithm that is based on an approximate 
series expansion. This algorithm was implemented in IBM's Power3™. Their algorithm provides 
no advantages over Goldschmidt's algorithm. In double precision, their algorithm requires 8 mul- 
tiplications and the longest chain of dependent multiplications consists of 4 multiplications. 

We present a version of Goldschmidt's division algorithm that uses directed roundings. We de- 
velop a simple general parametric analysis of tight error bounds for our version of Goldschmidt's 
division algorithm. Our analysis is parametric in the sense that it allows arbitrary one-sided errors 
in each intermediate computation and it allows an arbitrary number of iterations. In addition, we 
suggest four practical simplified settings in which errors in intermediate computations are not ar- 
bitrary. For each of these four settings, we present a closed error formula. The advantage of closed 
formulae is in simplifying the task of finding optimal parameter combinations in implementations 
of Goldschmidt's division method for a given target precision. 

We demonstrate the advantages of our error analysis by showing how it could lead to savings 
in cost and delay. For this purpose we consider Oberman's [28] floating-point micro-architecture 
used in the AMD-K7™ design. We present a micro-architecture that implements our version 
of Goldschmidt's algorithm and follows the micro-architecture described in [28]. The modules 
building our micro-architecture were made as similar as possible to the modules in [28]. This 
was done so that the issue of the precisions of the lookup table and multiplier could be isolated 
from other issues. Based on our analysis, we use a smaller multiplier (70 x 74 bits compared to 
76 x 76 in [28]) and we allow a slightly larger initial error (2" 13 - 51 compared to 2" 13 75 in [28]). 
Based on the cost models of Paul & Seidel [29] and Mueller & Paul [25], we estimate that our 
parameter choices for multiplier widths and initial approximation accuracy reduce the cost of the 
micro-architecture by 10.6% compared to the parameter choices in [28]. 

The paper is organized as follows. In Section C.2 we present Newton's method for division 
and then proceed by presenting Goldschmidt's algorithm as a variation of Newton's method. In 
Section C.3, a version of Goldschmidt's algorithm with imprecise intermediate computations is 
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presented as well as an error analysis. In Section C.4 we develop closed form error bounds for 
Goldschmidt's method with respect to four specific settings. In Section C.5 we present an alterna- 
tive micro-architecture to [28] and compare costs. 



In this section we present Newton's method for computing the reciprocal of a given number. We 
then continue by describing a version of Goldschmidt's division algorithm [16] that uses precise 
intermediate computations. We show how Goldschmidt's algorithm is derived from Newton's 
method. The error analysis of Newton's method is used to analyze the errors in Goldschmidt's 
algorithm. 

C.2.1 Newton's method. 

Newton's method can be applied to compute the reciprocal of a given number. To compute the 
reciprocal of B > 0, apply Newton's method to the function f(x) = B — 1/x. Note that: (a) the 
root of f(x) is which is the reciprocal we want to compute, and (b) the function f(x) has a 
derivative f'(x) = x~ 2 in the interval (0, oo). In particular, the derivative /'(x) is positive. 

Newton iterations are defined by the following recurrence: Let x 0 denote an initial estimate 
xq j=- 0 and define Xi+i by 



C.2 Goldschmidt's division algorithm 




= Xi • (2 - B • Xi) 



(C.l) 



A few iteration steps are visualized in Figure C.l. 
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Figure C. 1 : The first iteration steps of the Newton approximation of 1/5. 



Consider the relative error term defined by 



e i = ^ r - 



= 1 - B • Xi. 



It follows that 



= 1 — -B • 

= 1 - £ • x, • (2 - 5 • x,) 
= (1-B.* i ) 2 



(C.2) 



Equation C.2 has three implications: 

1. Convergence of x { to 1/B at a quadratic rate is guaranteed provided that the initial relative 
error is less than 1. Equivalently, convergence holds if x 0 e (0, f ). 



2. For i > 1, the relative error is non-negative, hence, x» < 1/5. This property is referred to 



114 



Algorithm 3 Goldschmidt-Divide(yl, B ) - Goldschmidt's iterative algorithm for computing A/B 

Require: |e 0 | < 1. 
1: Initialize: N- X <- A, £>_i <- B, F_i <~ 
2: for i = 0 to k do 

3: Ni «- Ni-x • 

4: A «- A-l * J^-l. 

5: F< 4- 2 - A- 

6: end for 

7: Return(A^) 

as onesided convergence, 

3. If B 6 [1, 2), then also the absolute error decreases at a quadratic rate. Hence, the number 
of "accurate" bits doubles in every iteration, and the number of iterations required to obtain 
p bits of accuracy is logarithmic in p. 

The disadvantage of Newton's iterations, with respect to a pipelined multiplier, is that each 
iteration consists of 2 dependent multiplications: a» = B • X* and ft = Xi • (2 — c*i). Namely, the 
product ft cannot be computed before the product a* is computed. 

C.2.2 Goldschmidt's Algorithm 

In this section we describe how Goldschmidt's algorithm can be derived from Newton's method. 
Here, our goal is to compute the quotient A/B. Goldschmidt's algorithm uses three values N u Di 
and Fi defined as follows: 

Ni = A-Xi 
Di = B-Xi 
Fi = 2 — A- 



Consider Newton's iteration: = Xi • (2 - B • Xi). We may rewrite an iteration by 
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which when multiplied by A and B y respectively, becomes 

D i+1 = Di Fi 

Since x { converges to l/B y it follows that AT* converges to A/B and £>< converges to 1. 

Note that: (a) Ni/Di = A/B, for every i; (b) Ni converges to A/B at the same rate that x t 
converges to 1/B; and (c) Let A,B>0. Since the relative error in Newton's is non-negative, 
for i > 1, it follows that TV, < A < 1 and F { > 1 for i > 1. 

As in Newton's iterations, the algorithm converges if |e 0 | < 1 and the relative error decreases 
quadratically. One could use a fixed initial approximation of the quotient. Usually a more accurate 
initial approximation of 1/B is computed by a lookup table or even a more elaborate functional 
unit (c.f. [6, 36]). 

Algorithm 3 lists Goldschmidt's division algorithm. Given A and B the algorithm computes 
the quotient A/B. The listing uses the same notation used above, and iterates k times. 

Observe that the two multiplications that take place in every iteration (in Lines 3-4) are inde- 
pendent, and therefore, Goldschmidt's division algorithm is more amenable to pipelined imple- 
mentations. Figure C.2 depicts and compares the flowcharts of the iterations of Newton's method 
and Goldschmidt's algorithm. The initial approximation is assumed to depend on the value of B. 
A closer look at Figure C.2 reveals that A; iterations of either Newton's method or Goldschmidt's 
algorithm require 2k + 1 multiplications (the unfilled boxes in the figure refer to initializations in 
which multiplication does not take place). These 2k -h 1 multiplication must be linearly ordered in 
Newton's method implying a critical path of 2k + 1 dependent multiplications. In Goldschmidt's 
algorithm the two multiplications that take place in every iterations are independent, hence the 
critical path consists only of k + 1 multiplications. 

An error analysis of Goldschmidt's algorithm with precise arithmetic is based on the following 
claim. 
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a) ^FpROX(l/B) 



Xo 



D 0 :=X 0 B (T) 



X\ := Xo • Fo 



*1 



X X 



Di :=Xt B (?) 
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N k :=X k A 
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D-i := B 




R. 1 :=APPROX(yB) 





D 0 ■= D-i ■ F-i 

Fg :=2- Dp 



Do 



1 



iteration -1 



N. 



\No := N-i ■ F-i 



£>x := D 0 • F 0 
ft ==a-X)i 
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N! := N 0 ■ F 0 
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Figure C.2: Sched5feof the Iterations of: a) Newton's method and b) Goldschmidt's division 
algorithm using an initial approximation for l/B. The numbers in circles indicate the sequence of 
the multiplications involved. A bound on the relative error e { of iteration i appears in each iteration. 



Claim 28 The following equalities hold for i > 0: 

A = 1 - eg 1 (C.3) 
Fi = 1 + cf (C.4) 

The key difficulty in analyzing the error in imprecise implementations of Goldschmidt's algorithm 
is due to the violation of the invariant Ni/Di = A/B. Consider the equality 

N k = A/B • £>o • F 0 • F x • ... • 

= A/S.(l-e 0 )-(l + e 0 ).(l + eg).(l + e 0 ). ... .(i + eg*- 1 )- 

Imprecise Z} 0 , F 0 , . . . , F fc _i accumulate to an imprecise approximation of A/B. 
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How many iterations are required? Consider an initial relative error of |e 0 | < 2" a , and assume 
that a quotient (reciprocal) with a relative error smaller than 2~ n is required. Implementations 
with precise intermediate computations of Newton's method and Goldschmidt's algorithm require 
k = \log 2 (n/a)'] iterations. 

Precision of intermediate multiplications. Let len(x) denote the number of bits in the binary 
representation of x. Precisions above 2n (where n denotes the number of bits used to represent each 
of the operands A and B) are usually considered too high for practical implementations. The fol- 
lowing claim deals with the lengths of the operands in intermediate calculations of Goldschmidt's 
algorithm. The claim shows that the required precision is (a + n) This renders implementations 
with precise arithmetic impractical. 

Claim 29 Let len(A) = len(i?) = n and len(F_i) = a. The lengths of the operands N u Di and 
Fi in Goldschmidt's algorithm are 2 i • (a + n), for i > 0. 

C.3 Imprecise Intermediate Computations 

This section contains the core contribution of our paper. We present a version of Goldschmidt's 
algorithm with imprecise intermediate computations. In this algorithm the invariant Ni/Di = A/B 
of Goldschmidt's algorithm with precise arithmetic does not hold anymore. We then develop a 
simple parametric analysis for error bounds in this algorithm. The error analysis is based on 
relative errors of intermediate computations. The setting is quite general and allows for different 
relative errors for each computation. 

We define the relative error as follows. 

Definition 7 The relative error ofx with respect to y is defined by 



y - x 



y 
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Note that one usually uses the negative definition (i.e., (x - y)/y). We prefer this definition since 
is helps clarify the direction of the directed roundings that we use. 

The analysis begins by using the exact values of all the relative errors. The values of the relative 
errors depend on the actual values of the inputs and on the hardware used for the intermediate 
computations. However, one can usually easily derive upper bounds on the absolute errors of 
each intermediate computation. E.g., such bounds on the absolute errors are simply derived from 
the precision of the multipliers. Our analysis continues by translating the upper bounds on the 
absolute errors to upper bounds on the relative errors. Hence we are able to analyze the accuracy 
of an implementation of the proposed algorithm based on upper bounds on the absolute errors. 

An interesting feature of our analysis is that directed roundings are used for all intermediate 
calculations. Surprisingly, directed roundings play a crucial role in this analysis and enable a 
simpler and tighter error analysis than round-to-nearest rounding (c.f. [28]). 

C.3.1 Goldschmidt's division algorithm using approximate arithmetic 

A listing of Goldschmidt's division algorithm using approximate arithmetic appears in Algo- 
rithm 4. The values corresponding to A^, and F { using the imprecise computations are denoted 
by jVjf, D\ and F{, respectively. 

Directed roundings are used for all intermediate calculations. For example, N- is obtained by 
rounding down the product N-_ x • F-_ v We denote by ra* the relative error of N{ with respect to 
Nl_ x • F(_ v Since N-_ x • F-^ is rounded down, we assume that ra* > 0. Similarly, rounding down 
is used for computing F{ (with the relative error /*) and rounding up is used for computing D\ 
(with the relative error d*). 

The initial approximation of the reciprocal 1/B is denoted by FL X . The relative error of FL\ 
with respect to 1/B is denoted by eo- We do not make any assumption about the sign of e 0 . 
Our error analysis is based on the following assumptions: 

1 . The operands are in the range A, B 6 [1, 2). 
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2. All the relative errors incurred by directed rounding are at most 1 /4. This assumption is 
easily met by multipliers with more than 4 bits of precision. 

3. We require that |eo| + 3do/2 + / 0 < 1/2. Again, this assumption is easily met if the multi- 
plications and the initial reciprocal approximation are precise enough. 

4. The initial approximation FL\ of 1/B is in the range [1/2, 1]. This assumption is easily met 
if lookup tables are used. 



Algorithm 4 Goldschmidt-Approx-Divide(^4, B) - Goldschmidt's division algorithm using ap- 
proximate arithmetic 

1: Initialize: N'_ x <- A 9 D'^ <- B y F'_ x <- ~ ~~ " 

2: for i — 0 to k do 

3: N ^(l^ ni) . NU . FU , 

5: - (1 - /<) • (2 - DQ. 

6: end for 

7: Return(iV t ') 



Definition 8 77*e relative error ofN- with respect to A/B is 

p{Ni) = A/B ' 

C.3.2 A simplifying assumption: strict directed roundings 

The following assumption about directed rounding used in Algorithm 4 helps simplify the analysis. 

Definition 9 (Strict Directed (SD) rounding) Rounding down is strict if x > 1 implies that 
rnd(x) > 1. Similarly, rounding up is strict if x < 1 implies that rnd(x) < 1. 

Observe that, in general, rounding down means that rnd(x) < x> for all x. Often the absolute 
error introduced by rounding is bounded by e > 0, namely x — e < rnd(x) < x. Strict rounding 
down requires that if x > 1, then rnd(x) > 1 no matter how close x is to 1. In non-redundant 
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binary representation strict rounding is easily implemented as follows. Strict rounding down can 
be implemented by truncating. Strict rounding up can be obtained by (i) an increment by a unit 
in each position below the rounding position and (ii) truncation of the bit string in positions less 
significant than the rounding position. 

Assumption 30 (SD rounding) All directed roundings used in Algorithm 4 are strict. 
C.3.3 Parametric Error Analysis 

Lemma 3 1 bounds the ranges of N- and F{ in Algorithm 4 under Assumption 30. This lemma is 
an extension of the properties Di < 1 and F { > 1 (for i > 1) of Algorithm 3. 

Goldschmidt already pointed out that since Fi tends to 1 from above, one could save hardware 
since the binary representation of Fi begins with the string 1.000 — An analogous remark holds 
for However, Lemma 31 refers to the inaccurate intermediate results (i.e., D\ and F{) rather 
than the precise intermediate results (i.e., Di and Fi). Parts 2-3 of lemma 31 show that the same 
hardware reduction strategy applies to Algorithm 4, even though intermediate calculations are 
imprecise. 

Definition 10 Define S it fori > 0, as follows: 

Si := < 



|e 0 | + 3d 0 /2 fori = 0 
Si-i + /i-i otherwise. 



Lemma 31 ( ranges of D\ and F[) 
The following bounds hold: 

1. D' 0 e [l-5 0 ,l + <*o] ^ (1/2,3/2). 

2. D\ e [1 - 5 iz 1], for every i > 1. 

3. F( G [1, 1 + Si], for every i > 1. 
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4. D' 4 < D' i+1 , for every i > 1. 

Proof: Part (1): Since D' 0 = (1 + d 0 ) ■ D'_ x ■ F'_ x , it follows that D' 0 = (1 + d 0 ) • (1 - eo) = 1 - 
e 0 + d 0 (l - e 0 ). Since d 0 > 0 and |e 0 | < 1/2, it follows that (1 - e 0 ) < 3/2 and 

I -S 0 < l-|e 0 | < < 1 + (|e 0 | + 3do/2) < 1 + S 0 . 

The bounds of 1/2 < l-<5 0 < 1+S 0 < 3/2 follow immediately from the condition 0 < 6 0 < |e 0 | + 
3d 0 /2 + /o < 1/2. 

Part (2): We prove that if A-i e [1 - 5^, 1 + 6^], then A € (1 - 1]. It then follows 
from Part (1) that £>i is in the required range, and then by induction, D* is in the required range, 
for every i > 1. 

Let U-i = A-i - 1- The assumption £>»-i e (1 - 1 + implies that \ti-i\ < 5i-i. 
The value of D\ can then be written as: 

D\ = (1 + • (1 + *<_!) • (1 - f^-x) • (1 - /<_!) 

= (i^-(i^^-^A_^ 

>1 <1 <1 

>(1 -*?_,) (! -/i-i) ( C5 > 

>(l-«?-i-/i-i) 

>(l-Ci-/i-i) 

= 1-* 

The multiplication by (1 + d*) in the second line of Equation C.5 models the strict rounding up of 
(1 - t}) ■ (1 - fi) < 1. Hence, D[ < 1, and Part (2) follows. 

Part (3): By Part (2), 1 - 6 t < D\ < 1, hence 1 < 2 - D' { < 1 + S t . Recall that F{ is obtained 
by strict rounding down of 2 — D\, hence 1 < F[ < 1 + Si. 

Part (4): From Part (3) it follows that F( > 1. Since D' i+l is obtained by rounding up £>< • F[, 



122 



it follows that D' i+l > D' iy as required. □ 
The following claim summarizes the relative error of Goldschmidt's division algorithm using 
approximate arithmetic. 

Theorem 32 For every i > 0, the relative error p{N!-) = A/ %^* satisfies 

m < p{N[) < tt, + 4, (C.6) 

where ir< is defined by 7r t - = 1 - (1 - n<) ■ > 0. 

Proof: We decompose p(N!) as follows. 

_ A/B - Nj _ A/B - Ag/CgU • FU) , jW-i • g-i) - rr _, 
P ™ A/B A/B + A/B K ' 

We first show that the first term on the right hand side equals 7r 4 . Observe that ^ can be expanded 
as follows: 



^A t~t 1 - 

J=0 J 



(C.8) 



Equation C.8 implies that N^/D^ is non-increasing. 



N' 

An expansion of D , ^ yields 



by definition of 



^-1 

t-1 1 

n ni 



Note that in particular ^ ^ ■ < ^, hence J is non-negative. It follows that 
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p{ jy ^} F , - ) satisfies 



o( N< = A/B-N'J{DU-FU) 
^DUF'J A/B 



1 — 71 j 



i-(i-*).n^ 



(C.9) 



= 7IY 



We now prove that the second term on the right hand side of Eq. C.7 is non-negative and 
bounded by Si. Consider the numerator 

N- N- 



(byeq.C.9) = (1 - w,) • | • (1 - • i^.J 
> 0. 

To complete the proof we only need to show that 0 < 1 — D'^ ■ F/_j < <5j. By Parts 1-2 of 
Claim 31, it follows that Z^_ x € [1 - 6^ u 1 + <5i_i] and Fl_ x = (1 - • (2 - D^). Hence 
• is bounded by 

l-<Ji < < D'i-iF'i-i < D' i _ l .{2-D' i _ l ) < 1, 

and the claim follows. □ 
For i = 0 it can be verified that \p(Nq)\ < ir Q + 5 0 . 

A somewhat looser (yet easier to evaluate) bound on the relative error follows from 

t i-l 

*<<x>i + ;C d *- (do) 

j=0 j=0 

The proof of equation CIO is as follows. Note that (1 + x)" 1 > 1 — x, for x e [0, 1) and 



> • 
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- «*) > 1 - Ei *i> for Xi € [0, 1). Hence, 

«-i-<i-«0-ICtT? 

< i - (i - m) . n^^ 1 - n ^ * ( i - d i) 

^ 1 -(i-(^+X^ i ^+*))) 

C.3,4 Deriving bounds on relative errors from absolute errors 

In this subsection we obtain bounds on the relative errors p{N , i ) from the absolute errors of the in- 
termediate computations. The reason for doing so is that in an implementation one is able to easily 
bound the absolute errors of intermediate computations; these follow directly from the precision 
of the operation, the rounding used (e.g., floor or ceiling), and the representation of the results 
(binary, carry-save, etc.). 

Consider the computation of A^. The relative error introduced by this computation is n u and 
N[ equals (1 - n*) • N-_ x • F{_ x . An accurate computation would produce the product N-_ x • F{_ x . 
Hence, the absolute error is • N<_ x • F{_ x . 

Definition 11 The absolute errors of intermediate computations are defined as follows: 

nepsi = rii • N^ x • F-_ x 
dep Si = di • • F[_ x 
fep Si ± /«-(2-D{). 

In an implementation, the exact absolute errors are unknown. Instead, we use upper bounds on 
the absolute errors. We denote these upper bounds as follows: neps i > nepsi, deps i > depsi and 
feps i > fepsi. 

The following claim shows how one can derive upper bounds on the relative errors from upper 
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bounds on the absolute errors. 

Claim 33 (from absolute to relative errors) If A, B G [1, 2), then for i > 2 the relative errors are 
bounded by: 

0 < rii < 2neps i /(l — iti-i — 
0 < di < depsJil-S^) 
0 < fi < f^Si 



Proof: It follows from Definition 1 1 that 

_ nepSi depsi fepsi 

By Theorem 32 it follows that, for i - 1 > 1, 

^-i>|-(l-iri-i-*-i)- 
By Part 3 of Lemma 31 it follows that, for i - 1 > 1, > 1. Hence 

rii < nepSi • — . 

1 — -Ki-i - di-i 

By Parts 2 and 3 of Lemma 31 it follows that 1 > D[_ r > 1 - S^i and F[_ x > 1. Hence the 
bounds on di and fi follow. □ 
A careful reader might be concerned by the fact that and 7T;_i appear in the above bounds 
on the relative errors and d^ When analyzing the errors, one computes upper bounds for all 
relative errors from the first iteration to the last. These bounds are used to compute upper bounds 
on Si-i and n^u which in turn are used to bound rii and d*. In section C.6, bounds are given for 
the relative errors in the first and second iteration (see Claim 37 in Section C.6). 
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C.4 Closed Form Error Bounds in Specific Settings 



In this section we describe four settings of the relative errors that enable us to derive closed form 
error bounds. The advantage of having closed form error bounds is that such bounds simplify the 
task of minimizing an objective function (modeling cost or delay) subject to the required precision. 
Closed form error bounds also enable one to easily evaluate the effect of design choices (e.g., initial 
error, precision of intermediate computations, and number of iterations) on the final error. 

C.4.1 Setting I: n*, di < n and fi = 0. 

Setting I deals with the situation that all the relative errors n i3 di are bounded by the same value n. 
In addition it is assumed in this setting that fi = 0, for every i. The justification for Setting I is 
that if all internal operands are represented by binary strings of equal length, then it is possible to 
bound all the relative errors n*, di by the same value. The relative errors fi can be assumed to be 
0, if the computations F( = (2 — D[) are precise. 

Using Theorem 32 and Equation C.10, the relative approximation error p(N!) = A/ ^~ iV ' in 
Setting I can be bounded by: 



C.4.2 Setting II: n u di<n and fi/5f is exponential in — k. 

In setting II it is assumed that fi/Sf < a • 2~ fc , where a is an appropriately chosen constant and k 
is the number of iterations. In addition, it is assumed that n*, < n, for every i. 

Claim 34 



0 < p(Nl) = 



A/B - 
A/B 



< (2< + l)-n + (|e 0 |+3n/2) 2 \ 



(l+a-2- fc ) : 



fc\2*-l 



<6?-e°. 
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Observe that in Setting II 5 k converges quadratically and that S k is larger than 5%" only by a 
constant factor which is independent of the initial approximation error. 

Based on Theorem 32 and Equation CIO the error bound in setting II satisfies: 

p(N' k ) < (2k + 1) • ft + e<* • (|e 0 | + 3rf 0 /2)) 2fc . 
C.4.3 Setting III: n u di < n and fi/5? is constant. 

In setting III it is assumed that fi/6f < C, where C is an appropriately chosen constant which is 
independent of the number of iterations k. In addition, it is assumed that n u di < h, for every i. 

Claim 35 

4 < (i + C) 2 *- 1 -^. 
Based on Theorem 32 and Equation C. 10 the error bound in setting in satisfies: 
p(N' k ) < (2k + 1) • h + ((1 + C) ■ (|e 0 | + 3d 0 /2)) 2 \ 

C.4.4 Setting IV: n u di<h and < f for every i. 

In setting IV the assumptions are: (i) n i5 d { < n, for every i, and (ii) < f < 1/8, for every i. 
Hence, Si < Sf_ x + / for alH > 0. 

The following claim bounds the error term 5 k corresponding to the fcth iteration of Algorithm 4. 



Claim 36 Let a = 




For every k > 0 the following holds: 



5 k < f+ m^a 2 ^ 1 - 2 -^, (a 2fc - 2 -5 0 2fc - 1 +/) 2 ,9/ 2 } 
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Proof: We choose the substitution Q = 5 t - /, for i > 0, and Co = <*o- Then, for i > 0: 

101 = \Si-f\ 

< £l 

= (G-1 + /) 2 (C.ll) 

Observe that if |C»— x| < "J~f, then \&\ < 2 • / < yff. This is proved as follows. Since / < 1/4, it 
follows that 2/ < y''/. By Eq. C.l 1 we get: 

ICil < (IG-1I + /) 2 

< (yff + f) 2 

< 2/. 

The last line holds since / < (\/2 - l) 2 « 0.171. 

Let t denote the highest index such that Coi • • • , O-i > y/J- Applying Equation C.l 1 we get 
for every i < £: 

ICil < (IO-1I + /) 2 

< (ICi-il + \/7-ICi-il) 2 

= (a-Ci-i) 2 - 



Hence by induction: 



101 < a 2 ' +l - 2 -C 0 2 ', (C.12) 



If I > k, then the claim holds. Otherwise, there are two possibilities: (i) £ < k — 2 or 
(ii) I = k — 1. If case (i) holds, then Q < ^//, and hence \Q\ < 2 • /, for every j > L Applying 
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Equation C. 1 1 we get: 

lai < (ia-ii+/) 2 

< 9/ 2 . 

If case (ii) holds then by applying Equation C.l 1 we get: 

I&I < (101 + /) 2 

< (a 2 *" 2 • 1 + /) 2 (byEq.C.12). 

Combining these cases (namely: £ > k,£ = k — 1 and £ < k — 1) yields: 

C fc < max{a 2t+l - 2 • 5* , (a 2 *" 2 • 1 + /) 2 , 9/ 2 }. 

and the claim follows. □ 
Note that a slightly looser bound that does not involve a max function can be written as: 

4 < / + a 2fc+1 - 2 A 2 * + 7/ 3/2 . 

Based on Theorem 32 and Equation CIO the error bound in setting IV satisfies: 

p(N' k ) < (2k + 1) • n + / + maoc{/ +1 - 2 . 8f , (a 2 ^ 2 • + f) 2 , 9/ 2 }. 

One can easily see that, due to the first term, there is a threshold above which increasing the 
number of iterations (while maintaining all other parameters fixed) increases the bound on the 
relative error. Moreover, the contribution of the error term / to p(N' k ) does not increase with the 
number of iterations k (as opposed to n). This implies that in a cost effective choice one would use 

f > h. 
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C.5 Application: An Alternative FP-DIV Micro-architecture 
for AMD-K7™ 

In this section we propose an alternative FP-DIV micro-architecture for the AMD-K7 micropro- 
cessor [28]. This alternative micro-architecture is a design that implements Algorithm 4. Our 
micro-architecture uses design choices that are similar to those of [28] to facilitate isolating the 
effect of precisions on cost. Our error analysis allows us to accurately determine the required 
multiplier precision and thus both save cost and reduce delay. 

Overview micro-architecture. The FP-DIV micro-architecture of the AMD-K7 microprocessor 
is described in [28]. The micro-architecture is based on Goldschmidt's algorithm. We briefly out- 
line this micro-architecture: (i) Round-to-nearest rounding is done in intermediate computations 
(as opposed to directed rounding suggested in Algorithm 4). (ii) The design contains a single 
76 x 76-bits multiplier. This means that the absolute errors neps* an d deps i are identical during 
all the iterations (i.e., since round-to-nearest is used, neps i = deps i — 2~ 76 ). However, our alter- 
native micro-architecture may use smaller multipliers (even multipliers in which the multiplicands 
do not have equal lengths) provided that the error analysis proves that the final result is accurate 
enough, (iii) Intermediate results are compressed and represented using non-redundant binary rep- 
resentation. This means that Assumption 30 on strict directed rounding is easy to implement in 
our alternative micro-architecture. (Recall that directed rounding is used in Algorithm 4.) (iv) The 
computation of F- is done using one's complement computation. This means that the absolute er- 
ror fepSi is identical during all the iterations, and that the error analysis of Setting IV is applicable 
for our alternative architecture, (v) Final rounding of the quotient is done by back multiplication. 
Our alternative micro-architecture uses the same final rounding simply by meeting the same error 
bounds needed in the final rounding of [28]. 
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Required final precisions. The micro-architecture in [28] supports multiple precisions: single 
precision (24, 8) in one iteration, double precision (53, 11) in two iterations, an extended precision 
(64, 15) and an internal extended precision (68, 18) in three iterations. Final rounding is based 
on back-multiplication: namely, comparing N' k ♦ B with A. In general, correct IEEE rounding 
based on back-multiplication requires that p(N' k ) < 2~( p+1 \ where p denotes the precision. (The 
description of the required precision for correct rounding in [28] is somewhat confusing since it is 
stated in terms of a two sided absolute error. For example, the absolute error in the 68-bit precision 
is bounded by 2~ 70 .) 

To summarize, the upper bounds on the relative errors are as follows: (i) for single preci- 
sion: p(N{) < 2~ 25 , (ii) for double precision: p(N^) < 2" 54 , (iii) for extended double precision: 
p(N£ < 2" 65 , and (iv) for the 68-bit precision: p(N£) < 2" 69 . 

Note that the bound for the 64-bit precision is weaker than the bound for the 68-bit precision. 
The bound for single precision is easily satisfied by the parameter choices needed to satisfy the 53- 
bit precision. Hence we focus below on two iterations for double precision and on three iterations 
for the 68-bit precision. 

From relative errors to multiplier dimensions. We first discuss the lengths of the multiplicands 
in Algorithm 4 since these lengths determine the cost and delay of the multiplier. The lengths (or 
precisions) of the multiplicands is derived from upper bounds on the relative errors n i5 di and 

In Algorithm 4 the multiplier performs the multiplications: N- • F[ and jD< • F(. It is reasonable 
to assume that one multiplicand is used for N[ and D\ and that the other multiplicand is used for 
F[. (We later elaborate on which multiplicand should be Booth recoded.) 

For simplicity, let us assume that Setting IV holds (i.e., n*, d» < ft and fa < /, for every i). 

We start by deriving the required length of the first multiplicand (used for N[ and D^). Let a 
denote the length of the first multiplicand; this means that this multiplicand is represented by a 
binary string with bits in positions [0 : (a - 1)]. Since the product is rounded (either down or up) 
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to length a, the absolute errors satisfy: 



\nepsil \depsi\ < 2~ (a ~ 1) . 



Let us first focus on the constraint < n. By Claim 33 it follows that for i > 2: 



0 < Hi < 



42" a 



1 — 7T t _i — Si-i 



1 — TTi-i — 5,-1 



Hence it suffices for a to satisfy: 



4-2" a 



< fi. 



1 — 7Ti_i — 



Which implies 




1 - 7Tj_x - 



1 



Similar lower bounds for a are derived for the cases i = 0, 1 using Claim 37. Hence a should be 
slightly larger than log 2 (l/n) + 2. 

Applying Claim 33 with respect to di yields that D[ needs to be represented using slightly more 
than 1 + log 2 (l/n) bits. Hence, supporting the relative error requires an extra bit compared to 
d t . 

We remark that this extra bit stems from the fact that N[ can be less than 1. One could apply 
normalization shifts to the binary string that represents N- to reduce, in effect, the absolute error 
when N- < 1. This implies a reduction of the length a of the first multiplicand by 1 at the cost 
of adding the circuitry required for the normalization shift. Moreover, the normalization shift can 
take place in the 2nd pipeline stage (which has some slack with respect to delay) and help reduce 
the delay of the 3rd pipeline stage (which contains the addition tree of the multiplication tree). 
We do not take this optimization method into account in our comparison since it does not follow 
directly from our error analysis (which is the main contribution of this paper). 

Applying Claim 33 and Claim 37 with respect to /» yields (due to i = 0) that the length of the 
second multiplicand should be greater than or equal to log 2 (l/ f) + 1 + log 2 (l/(l — 6 0 )). 
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Optimizing the error parameters. Given a relative error bound on p(N' k ) 9 a combination of 
relative errors for intermediate computations is feasible if our error analysis shows that the relative 
error p(N' k ) is smaller than the required bound. In this paragraph, we search for feasible combi- 
nations that minimize the sizes of the multiplier and lookup table. We used a cost function that is 
based on the cost model of Paul & Seidel [29] for Booth Radix-8 multipliers and the cost model 
of Paul & Mueller [25] for lookup tables, adders, etc. (Formulas for hardware costs appear in 
Section C.7.) 

We demonstrate the power of our error analysis in finding optimal error parameters using Set- 
ting IV. Three error parameters determine the relative error p(N f k ) in Setting IV: n, f and e 0 . 
Figure C.4 depicts the results of a three dimensional search for the triple (n, /, e 0 ) that leads to the 
cheapest design that supports IEEE double precision (53-bits). Feasible pairs (n, /) that support 
double precision division are depicted in the colored region above the curve in Figure C.4. Each 
point in this region is assigned a maximum value of e 0 , so that the triple (n, /, e 0 ) is a feasible 
combination of error parameters. The values of n and / determine the multiplier dimensions (and 
hence its cost), and the value of e 0 determines the size of the lookup table (and hence its cost). In 
this fashion we are able to associate a cost to every feasible pair (n, /) (we used here a smooth 
continuous version of the cost function so that tradeoffs are easier to interpret). Closed curves 
are used to connect parameter combinations that lead to designs with equal cost. The parameter 
combination corresponding to the cheapest double precision design is depicted in Figure C.4. The 
optimal parameter combination is - log 2 (/) = 55.67, - log 2 (n) = 57.74, and - log 2 (e 0 ) = 13.92. 

To simplify the search for the 68-bit internal precision we proceeded as follows. Setting I was 
employed to compute a good choice for e 0 . Setting I assumes not only a bound n for every rii 
and di\ it also assumes that = 0. We note, however, that following [28] our micro-architecture 
uses one's complement subtraction, so fi^O 1 . Nevertheless, Setting I is useful for the purpose of 
evaluating e 0 . Using Setting I, we compute feasible pairs (n, e 0 ), so that the relative errors satisfy 

! One may use injections in the micro-architecture to account for the missing ULP due to one's complement sub- 
traction. However, we do not pursue this correction because we wish to stick to a micro-architecture as similar as 
possible to that in [28]. 
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p(N^) < 2 -68 and p(N^) < 2" 53 . Figure C.3 depicts the curves of feasible pairs (n, e 0 ) for double 
precision and internal 68-bit precision, respectively. The area above these curves constitutes the 
feasible area per constraint, hence we are interested in the intersection of the areas above the curves. 
We conclude that the pair h = 2~ 70 81 e 0 = 2" 13 51 is a feasible pair for all the required precisions. 
It is also evident that it is of no use to try to increase e 0 above 2" 13 51 . (We do need to decrease h 
due to the fact that / ^ 0). 

Based on the analysis above, we fix e 0 = 2~ 13 51 , and search for an optimal feasible pair (n, /) 
for the 68-bit precision. Such a feasible pair is also a feasible pair for the remaining precisions. 
The top figure in Figure C.5 depicts feasible pairs (n, /) for 68-bit precision. Recall that the length 
of the first multiplicand length is slightly more than 2+log 2 (l/n), and that the length of the second 
multiplicand is slightly more than 1 + log 2 (l//)- One can show that this small quantity that must 
be added when determining the multiplicands' lengths is less than 0.09. The bottom figure depicts 
the cost of the micro-architecture as a function of / (when h is taken to be as large as possible). 
The cost function took into account the discretization due to computing the multiplicands' lengths. 
Multiplier dimensions are depicted in the bottom figure. Note that either multiplicand could be 
Booth recoded (we use radix-8 as in [28J). Each instance of multiplier dimensions (e.g., 70 x 74) 
means that the first multiplicand A is preprocessed to compute 3A, and the second multiplicand is 
Booth recoded. This explains the small effect that swapping the multiplicands has on the cost (e.g., 
70 x 74 vs. 74 x 70). The minimum cost is obtained for the feasible pair — log 2 n = 71.91 and 
— l°g2 / — 68.9. We conclude that multiplier dimensions 70 x 74 combined with a relative error 
bound e 0 < 2~ 13 51 are a feasible choice of error parameters 2 . These parameters lead to a savings 
in cost of 10.6% compared to the micro-architecture described in [28]. 

2 If one considers delay (rather than cost), then multiplier dimensions 74 x 70 lead to a slightly faster design 
(although slightly more expensive). Reduction in delay is due to the number of Booth digits: 24 Booth digits (corre- 
sponding to a 70-bit multiplicand) vs. 25 Booth digits (corresponding to a 74-bit multiplicand). Note that in [28] there 
are 26 Booth digits. 
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Figure C.3: Feasible parameter combinations of eo and h for double precision and extended double 
precision based on Setting I. 
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Figure C.4: Feasible (n, /) pairs using Setting IV for double precision are located in the shaded 
region above the curve. The closed curves depict pairs that lead to designs with equal costs. The 
central point depicts a feasible pair that leads to the cheapest design. 
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Figure C.5: Top: Feasible (n, /) pairs using setting IV for 68-bit precision (when e 0 = 2 -13 51 ). 
Bottom: Cost of design as a function of /. 
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C.6 Bounds on relative error for first two iterations 

Claim 37 

_ nepso ^ 2 ■ neps 0 
no — 7- — <. 



d 0 = 
fo< 



(l-eo) A/B ~ l-e 0 
depso 

(1 " eo) 
fepsp 

(I -So) 



ni < nepsi < nep Sl 

di < 



A/B ■ (1 - 7T 0 - 6 0 ) ■ i - (1 - /„) - (1 - /o) • (1 - 7T 0 - <5 0 ) 



(l-/o)(l-^) 
A < /epsi 

C.7 Cost model equations 

In this section we present the formulas that are used in the cost model. These formulas are on 
the cost models described in [25, 29]. Implementation of recoding and selection logic for Booth 
recoding radix-8 are from [4]. The cost formulas are listed below: 

cCLA{n) := 24 • fn] - 12 • |7o 52 (M)l + 6 
caddertree(n, t, S) := (n - S) • (t - 2) + (\log 2 (t)] - 2) • t ■ 8/2 
cBooth8tree(n,m) := 14 • addertree(\n + 8], 

3 

cBoothSrec(n, m) := 22 ■ \ m+ 

3 

cBooth8sel(n,m) := 18 - fn + 8] - f m + ^ 

3 

cBoothS(n, m) := cCLA(n + 2) + cBooth8tree(n, m) + cBooth8rec(n, m)+ 
+ cBooth8sel(n, m) + cCLA(n + m) 
Crom{A, d) := 0.25 • (A + 3)(d + log 2 (log 2 (d))) 
TLUibit) := Crom(2 Wt , bit - 1) 
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Appendix D 

Deeply pipelined multiplicative division 
with IEEE rounding using a full size 
multiplier with redundant feedback 

(A preliminary version of this paper appeared in [12]) 

Abstract 

We present a pipelined implementation of Goldschmidt's division algorithm with IEEE rounding. 
The core of the design is a Booth radix-8 multiplier, the operands of which are slightly longer 
than the target precision (i.e., a "full-size" multiplier). The pipeline has four stages: (i) initial 
approximation, (ii) computation of the 3x multiple of one operand and Booth recoding of the 
second operand, (iii) addition of the partial products to derive a carry-save representation of the 
product, and (iv) 2:l-addition and rounding selection. 

If the precision of the initial is roughly p/2 y where p denotes the target precision, then the 
latency is 9 cycles and a new division can be started after 6 cycles. This setting is especially suited 
for single precision (i.e., p = 24). In the case of double precision (i.e., p d = 53), if the precision 
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of the initial is roughly p d /4, the latency is 11 cycles and the restart-time is 8 cycles. 

Our implementation uses several techniques to reduce the latency, among them: (i) Both 
operands of the Booth multiplier can be either represented as non-redundant binary-numbers or 
redundant carry-save numbers, (ii) A novel rounding procedure, called dewpoint rounding, is in- 
troduced. Dewpoint rounding simplifies the rounding decision, unifies rounding for all rounding 
modes, and avoids using tables, (iii) Injection based rounding is used to implement directed round- 
ing of intermediate results. 

The paper is self-contained and includes a detailed and complete description for single preci- 
sion division. A simple error analysis is presented from which the exact multiplier dimensions and 
precision of the initial approximation are derived. 

D. 1 Introduction 

Goldschmidt's division algorithm was developed in the early 60's to overcome the problem of two 
dependent multiplications per iteration in Newton's method. The advantage of two independent 
multiplications per iteration is that these two multiplications can be parallelized or pipelined to 
reduce latency. Although Goldschmidt's algorithm requires fewer cycles compared to Newton's 
method, very few designs are based on Goldschmidt's algorithm. The main reason that designers 
avoided Goldschmidt's algorithm is that it is not self-correcting. Namely, inaccuracies in inter- 
mediate computations accumulate and cause the computed result to drift away from the accurate 
quotient. 

Recently, a general and parametric error analysis of Goldschmidt's algorithm was pre- 
sented [13]. This error analysis allows different one-sided errors in each iteration and deals with an 
arbitrary number of iterations. In [13], the implementation of Goldschmidt's algorithm in AMD's 
K7 floating-point microarchitecture [28] was revisited using this parametric error analysis. It was 
shown that AMD's K7 floating-point microarchitecture is quite conservative. Namely, the same 
IEEE rounded quotient can be computed with a smaller multiplier (i.e., shorter operands) and less 
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precise initial approximation (i.e., smaller tables). 

In this paper we address the issue of further improving the floating-point microarchitectures 
that implement Goldschmidt's algorithm. For a target precision of p bits, we focus on a microar- 
chitecture with the following characteristics: 

1 . The core of the microarchitecture is a "full-size" multiplier, namely, the dimensions of the 
multiplier are (p+£i) x (p+£2). The multiplier is a Booth radix-8 that is split into 3 pipeline 
stages. Such a pipeline allows for short clock cycles. 

2. The precision of the initial approximation is roughly p/2 bits. Since the number of bits 
of precision roughly doubles per iteration of Goldschmidt's algorithm, only one iteration is 
needed to meet the target precision. 

The exact multiplier dimensions and precision of the initial approximation are derived and min- 
imized using a simple error analysis provided in this paper. We present an implementation of 
Goldschmidt's algorithm with such a microarchitecture that has a latency of 9 cycles and enables 
starting a new division operation after 6 cycles. 

In [12] a microarchitecture that integrates a single precision and a double precision floating- 
point divider is outlined. To reduce cost and the clock period, the core of the design in [12] for 
double precision division is a half-size multiplier (i.e., somewhat larger than a p d + ei x (pa/2) + e 2 
multiplier, where pd = 53). Two iterations of Goldschmidt's algorithm suffice for double pre- 
cision since the precision of the initial approximation is roughly Such a microarchitecture 
constitutes, with respect to single precision (i.e., p s — 24), a microarchitecture with a full-size 
multiplier and the precision of the initial approximation is at least p 5 /2. Hence, this paper provides 
a complete and self-contained description of the microarchitecture for the case of single precision. 

Related work. Many alternative approaches have been considered for floating-point division, 
some of these approaches are based on non-multiplicative techniques [27, 24]. Comparisons be- 
tween multiplicative and non-multiplicative techniques are rare. Such a comparison appears, how- 
ever, in [24], where various non-multiplicative methods are compared with a multiplicative imple- 
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design 


latency (restart time) 
single precision 


latencv ( restart timel 
double precision 


JIM, ^UuoH/ loUlA 




o on* C?R(W 


SRT, (basic radix-4) 


78* (78) 


170* (170) 


SRT, (basic radix- 16) 


56* (56) 


120* (120) 


SRT, (short reciprocal) 


65* 


90* 


based on fpmult 


104 (85) 


160(136) 


VHR (very high radix) 


44* 


65* 


Boost VHR 


44* 


65* 


Proposed Implementation of Goldschmidt's Algorithm 


45 (30) 


66 (48) 



Table D.l: Latency comparison on division algorithms for single and double precision. The delay 
unit is the delay of a full adder. The top rows follow the model and estimates from [24]; the bottom 
row is the latency of the the proposed design. Latency estimates marked above with * require an 
additional delay for the IEEE rounding decision and selection and the implementation of all four 
IEEE rounding modes. 



mentation based on Newton-Raphson's method. Delays are expressed in [24] in terms of the delay 
of a full adder, and the timing model does not consider routing delays. However, these timing 
estimates enable a technology independent comparison. We extend the timing estimates provided 
by [24] for double precision to single precision as well. These timing estimates appear in Table D. 1 
together with an estimate of the delay of the implementation of Goldschmidt's algorithm presented 
in this paper. 

The estimate of the delay of our implementation is based on assuming that the clock period 
equals the delay of 5 full-adders in single precision and 6 full-adders in double precision. Since 
the latency is 9 and 11 cycles for single and double precision, respectively, the total delay is 45 
and 66 full-adders, respectively. Although this timing analysis is very rough (since it ignores wire 
delays and gate sizing), it indicates that the proposed implementation of Goldschmidt's division 
algorithm is rather competitive with the best floating point division algorithms developed to date. 

Techniques. The analysis presented in [13] assumes that all intermediate products are repre- 
sented in a non-redundant representation. To reduce latency, our microarchitecture permits carry- 
save representation of the multiplication operands. We therefore begin with a simple and short 
error analysis for one iteration of Goldschmidt's algorithm (see Sec. D.3). 
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We present a novel rounding procedure for IEEE floating point division (see Sec. D.4). We 
refer to this rounding procedure as dewpoint rounding. The procedure relies on an error range of 
the quotient that allows for only two candidates for the final IEEE rounded result. We associate 
with each candidate r a rounding interval I r . The rounding interval I r is simply the set of numbers 
that are rounded to r. The number that separates the two rounding intervals is called the dewpoint. 
The rounding decision is obtained by comparing the dewpoint and the exact quotient by applying 
back-multiplication. We present a unified dewpoint rounding procedure for all rounding modes 
and avoid rounding tables. Previous approaches for IEEE rounding in division compute a rounding 
representative of the exact quotient (i.e., a number that belongs to the same rounding interval that 
the exact quotient belongs to) and the round the representative [34]. 

We employ additional techniques, among them: (i) A Booth recoded multiplier that can be 
fed by either non-redundant binary operands or by redundant carry-save operands. This technique 
when applied to a Booth radix-8 multiplier enables reducing the feedback latency to two cycles. 
Previously, Booth multipliers with one operand in redundant carry-save representation were stud- 
ied [40, 8, 35]. (ii) Injection based rounding is used to implement directed rounding of intermediate 
results [9]. (iii) The rounding decision is based on the sign of the difference between the dewpoint 
and the exact quotient. We simplify the back-multiplication and the comparison by utilizing the 
fact that the difference between the dewpoint and the exact quotient has a small absolute value. 

Organization. In Section D.2, notation is introduced and the version of Goldschmidt's algorithm 
from [13] is reviewed. In Section D.3, an error analysis of this division algorithm appears. In 
Section D.4, dewpoint rounding is presented. In Section D.5, we list various techniques employed 
to save cost and delay in a hardware implementation. In Section D.6, we present an implementation 
of the division algorithm that is based on a full-size multiplier. For the sake of concreteness, the 
presentation is for single precision division. 
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D.2 Preliminaries 

Notation Let XiX i+ i • • - Xj € {0, 1}* denote a binary string. We often denote this string also by 
x[i : j]. We also refer to as x[i]. Since we deal with fractions (mostly in the binade [1, 2)), the 
weight associated with the bit x { is 2~\ Namely, a fraction is represented by the binary digit string 

X0.X1X2 - - ■ x p -i. 

Let a = (Ti(Ti+i — Gj G {0, 1, 2}* denote a carry-save encoded digit string (i.e. 0* € {0, 1, 2}). 
A carry-save encoded digit string a[i : j] is represented by two binary vectors <x'[i : j] and a"[i : j]. 
Each carry-save digit a[£] satisfies a[£] = a'[£] + &"[(]. 

Given a binary string x[i, j] and a carry-save encoded digit string a[i : j], we denote the values 
represented by these string by \x[i : j] \ and \a[i : respectively. 

Consider three binary vectors x y y,z. We denote 3:2-addition by FA(x,y,z). Namely, 
FA{x, y, z) means that the three binary vectors x, y, z are fed to a line of full-adders. The out- 
put of FA(x, y, z) is two binary vectors s and c (i.e., a carry-save number) that satisfy \x\ + \y \ + 

N = N + |c|. 

By truncating a carry-save encoded digit string a[i : j] after position q we simply mean chop- 
ping off the tail and leaving only a[i : q]. We denote truncation of a after position q by \cr\ q . 

We often regard a carry-save encoded digit string a also as two binary vectors a 1 and a". 
Hence if a is a carry-save encoded digit string and x is a binary vector, then FA(a, x) simply 
means FA(a\ a", x). 

Division Operation We consider the following task which captures the main difficulty in IEEE 
floating-point division. The dividend and the divisor are represented by p-bit binary strings A[§ : 
p— 1] andi?[0 : p— 1], where \A\, \B\ G [1,2). Our goal is to compute the binary string Q[0 : p — 1] 
that is the binary encoding of the rounded value of the normalized quotient. More formally, 

1. Let A = A if \A\ > \B\, and A = leftshift(A) if \A\ < \B\. This way, the quotient \A'\f\B\ 
is normalized (i.e., in the binade [1, 2)). 

2. h*q=\A'\f\B\. 



145 



3. Round q (according to the IEEE standard). The binary string Q[0 : p - 1] should satisfy 
|Q| - rnd(q). 

In this paper we consider single precision, namely p = 24. Four rounding modes are defined 
in the IEEE standard: round towards zero round towards zero (RZ), round to nearest even (RNE), 
round towards infinity (RI) and round towards minus infinity(RMI). 

Goldschmidt's division algorithm In [ 1 3] a variation of Goldschmidt's division algorithm based 
on directed rounding is presented and analyzed. The algorithm is listed below as Algorithm 5. 



Algorithm 5 Goldschmidt-Approx-Divide(A / , B) 
proximate arithmetic 

1: Initialize: NL X <- A\ D'_ x <- B, FL X <- 
2: for i = 0 to k do 

4: D ^ {l+di) , D[i . Fl ^ 
5: ^^(l-/ i ).(2- J D0. 

6: end for 

7: Return(JV;) 



- Goldschmidt's division algorithm using ap- 



Directed roundings are used for all intermediate calculations. For example, N- is obtained by 
rounding down the product N-^ F(_ v We denote by n» the relative error incurred when JV/_ 2 • F(_ x 
is rounded down. Rounding down translates to the assumption that n, > 0. Similarly, rounding 
down is used for computing F{ and rounding up is used for computing D\, 

In this paper we focus on an implementation of Algorithm 5 with k = 1 iterations. Namely, N{ 
is used as the estimate of the exact quotient for computing the rounded quotient. 



D.3 Error analysis 

In this section we prove bounds on the relative error of N[. 
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Definition 12 The relative error ofN[ with respect to ^| is denoted by p(N[) 

p W ~ \AV\B\ 

The relative errors n 0 , ni , d 0 , /o, e 0 are not known to the Algorithm 5. Instead, we are able to derive 
upper bounds on these relative errors (see Sec. D.6.3). Namely, n 0 > n 0i hi > n ly do > d 0 , fo > 
/o, and e 0 > |e 0 |. 

The following claim summarizes the bounds on p(N[). Note that all the assumptions on the 
relative errors are very easy to satisfy. 

Claim 38 Let 8 0 = |e 0 | + § • do- // 

0 < n 0 < n 0 < 1 0<ni<n!<l 0 < d Q < d 0 < 1 

0 < f 0 < fo < 1 0 < |eo| < eo < 1/2 6 0 < 1, 



then 



Proof: 



0 < p(N[) < 1 - (1 - n x ) • (1 - no) ■ (1 - /o) • (1 - eg - 4> ■ (1 + e 0 ) 2 ) 



Wj=(l-n 1 )-Ab-F 0 

= (1 - m) • (1 - n 0 ) • - ij^p ■ (1 - /o) • (1 + e 0 - d 0 + e 0 d 0 ). (D.l) 

To prove the lower bound we need to prove that N[ < Since 

0 < (I -n x ) • (1 -n Q ) • (1 - f 0 ) <h 



it suffices to show that 

0 < (1 - e 0 ) • (1 + e 0 - do + e 0 d 0 ) < 1- 



(D.2) 
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The assumption that |e 0 | < e 0 < 1/2 implies that e 0 - d 0 + e 0 do > — e 0 — |d 0 = —So- The 
assumption on S 0 and e 0 implies that the lower bound in Eq. D.2 follows. The upper bound in 
Eq. D.2 follows from the assumption that d 0 > 0 since 

(1 - e 0 ) • (1 + e 0 - d 0 + e 0 d 0 ) = 1 - (eg + <fc • (1 - e 0 ) 2 ) < 1. (D.3) 

We now prove the upper bound on p(N[). The assumptions imply that (1 — n 0 ) > (1 — ra 0 ) > 0, 
and the same holds for n\ and / 0 . By Eq. D.l, it suffices to prove that 

(1 - eo) • (1 + eo - do + e 0 d 0 ) > 1 - eg - d 0 • (1 + e 0 ) 2 > 0. (D.4) 
Equation D.4 follows from Eq. D.3 and the assumptions. The upper bound follows. □ 

D.4 EEEE rounding of the quotient 

In this section we describe how IEEE rounding is done in our algorithm. We present a novel 
method called dewpoint rounding. 

D.4.1 Background on IEEE rounding 

The IEEE Floating-Point Standard 754 [17] requires the implementation of all four IEEE rounding 
modes for all basic arithmetic operations and for all precisions supported. The four IEEE rounding 
modes comprise of round towards zero (RZ), round to nearest even (RNE), round towards infinity 
(RI) and round towards minus infinity(RMI). 

Back-multiplication is the common method for computing the IEEE rounded quotient. In this 
method, the approximate quotient is multiplied by the divisor B 9 and then this product is compared 
with the dividend A' [19, 34]. The rounded result is then computed to be within one ULP of the 
approximate quotient depending on the rounding mode, the sign, and the comparison. 
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We follow the method of reducing the four IEEE rounding modes to three rounding modes RZ, 
RI, and RNU (round to nearest up) based on the numbers' sign [32]. 

Observation 39 The exact quotient |j4'|/|i?| cannot be a midpoint between two representable 
numbers. Therefore RNU and RNE are equivalent rounding modes with respect to division. 

Based on the above reductions and Observation 39, we only need to consider the three rounding 
modes RZ, RI, and RNU. 

Injection based rounding [9] was introduced to further simplify rounding. Injection based 
rounding reduces these three rounding modes to truncation. Loosely speaking, this reduction is 
obtained by adding an injection that only depends on the rounding mode. Dewpoint rounding also 
uses injections, as described in the next section. 

D.4.2 Dewpoint Rounding 

Overview 

The concept of dewpoint rounding is inspired by the following concept of a dewpoint in meteorol- 
ogy: the dewpoint is the threshold temperature at which condensation of water in air occurs. 

If the quotient's estimate is close enough to the accurate un-rounded quotient, then rounding 
only needs to select between two values. We denote by RCi < RC 2 the two candidate values for 
the IEEE rounded result. Note that RC 2 = RCi + 2 _p+1 . For each rounding candidate RCj, we 
denote by Ij the rounding interval that comprises the pre-image of RCj with respect to rounding. 
Namely, Ij is the set of numbers that are rounded to RCj. The definition of RZ, RNU, and RI 
rounding implies that since RCi and RC 2 are successive representable numbers, the rounding 
intervals I\ and I 2 share an endpoint. This common endpoint is called the dewpoint. The rounded 
quotient is selected to be either RCi or RC 2 based on the sign of \A'\ — \B\ ■ dewpoint. 

Dewpoint rounding requires that there are only two rounding candidates. We are able to meet 
this requirement if the following assumption on N' k holds. 

Assumption 40 The quotient's estimate N' k satisfies: 0 < p(N' k ) < 2~ p . 
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Observe that our assumption is weaker than the conventional requirement of 0 < p(N' k ) < 2 p 2 
Computation of rounding candidates and dewpoint 

Given a carry-save encoding 9^[0 : t] of N' k , we wish to compute the dewpoint and the rounding 
candidates RC\ and RC^. A unified computation method based on injections is used for all round- 
ing modes. We begin by dealing with RZ rounding and then reduce the two other rounding modes 
RNU and RI to RZ. 

For the sake of simplicity we refer to the binary representation of N k and denote it by N' k [Q : t). 
Namely, N k = \N k [0 : t]\ = \%[0 : t]\. Note that full compression of N k is not required in the 
implementation. We use the binary representation just to simplify the description of the rounding. 
We use some notation related to the target precision p. 

Definition 13 We refer to multiples of2~ p + l as representable significands and to odd multiples of 
2~ p as mid-points. 

RZ rounding. Since 0 < p(N' k ) < 2~ p and ffl < 2, it follows that e [N' k , N' k + 2' p+1 ). Let 
LJVfcJp-i = \N' k [0 : p — 1]|. Truncation of N' k [Q : t] after position p - 1 yields: 

^e[L^ 1> L^J P - 1 +2-^ 1 ). 

We conclude that 

||f € [L^Jp-i, L^Jp-i + 2 • 2-*«). (D.5) 

In the RZ rounding mode the interval [[NjUp-i, L^Jp-i + 2 * 2 " p+1 ) is the union of exactly two 
rounding intervals: [L^Jp-i, IXJp-i+2" p+1 ) and [LWJUp-i+2 -1 * 1 , [N' k \ p - 1 +2-2- p + 1 ). Hence, 
the dewpoint in this case is [N k \ p -i + 2~ p+1 . The rounding candidates are simply RCi = I^/Up-i 

and i?C 2 = L^fcJ P -i + 2" p+1 - 

This method fails for RI and RNU since in these rounding modes the interval 
[l N k\p-ii l N h\p-i + 2 • 2~ p+1 ) intersects three rounding intervals. 
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Reduction of RI and RNU to RZ. We reduce the rounding modes RI and RNU to RZ by using 
injections [11] as follows. 

Let in J rnd (mode) denote an injection constant that depends only on the rounding mode, the 
target precision p, and the precision t of N' k . We define IN J md(mode) as follows. 



lNJ r „d(raode) = < 



0 if mode — RZ 

2~ p if mode = RNU (D.6) 

2-P+1-2"* if mode = RI. 



For x G [1, 2), the addition of IN J rnd (mode) reduces RI and RNU to RZ in the following sense 
(c.f.[ll]): 

RNU(x) = RZ(x + lN3 rnd (RNU)) 

(D.7) 

RI(x) = RZ(x + \NJ rnd (RI)). 



Hence it suffices to compute RZ(^ + iNJmd(mode)). 



RI & RNU Rounding. We now repeat the analysis for rounding modes RNU and RI as follows. 
Define 7V INJ as follows: 

~ N' k + lNJ rnd (mode). 

The following claim shows that the exact quotient is contained in one of two rounding intervals. 
We leave the proof as an easy exercise. 

Claim 41 



J L c / 

\B\ 



[L^.n,J p -i - 2-p L^.niJp-i + 3 • 2-f) if mode = RNU 
(L^.niJp-x - L^V INJ J P -i + 2-P+ 1 ) if mode = RI 



Note that for rounding mode RZ we have iV, NJ = N' k and the previous analysis for RZ is also 
covered by our current formalism. 
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The following claim summarizes dewpoint rounding for all three rounding modes. 



Claim 42 Let 



rnd mode 



\\B\J 



Rd = WnJp-i 

RC 2 ± + 

Rd + 2" p+1 if mode = RZ 
dewpoint = { Rd + 2~ p if mode = RNU 
RC\ if mode = RL 

For the rounding modes mode E {RZ,RNU,RI}, the IEEE rounded result has the following prop- 
erty: 

RC i i f(^ < dewpoint) 

RC i if ( ^ = dewpoint AND mode = RI), ( D - 8 ) 
HC2 otherwise. 

Note that the decision does not require division since < dewpoint) if and only if \A'\ < 
dewpoint • \B\. 

In the description of the division algorithm a unified notation is used. For this purpose we 
introduce the dewpoint displacement constant Cdew(™>ode) that is defined as as follows. 

2-( p " 1 > if mode =RZ 
Cde W (mode) = { 2" p if mode =RNU 
0 if mode =RI. 



The dewpoint satisfies: 



dewpoint = L^injJp-i 4- Cdswimode). 



152 



D.5 Hardware Optimizations 

D.5.1 Redundant Feedback & Partial Compression 

Addition trees in multipliers are not amenable to pipelining. Short clock cycles are therefore not 
achievable with reasonable cost if the addition tree has too many rows. Booth radix-8 recoding 
reduces the number of rows in an addition tree from n to l" 1 ^! . 

Booth radix-8 multipliers are usually implemented using a 3-stage pipeline, as follows: (i) 
precompute the 3x multiple of the first operand of the multiplier and recode the second operand, 
(ii) addition tree that computes a carry-save representation of the product, and (iii) final carry- 
propagate addition. 

Goldschmidt's algorithm performs only 2 multiplications per iteration. Hence running Gold- 
schmidt's algorithm on a 3-stage pipeline creates un-utilized cycles (i.e. "bubbles") in the pipeline. 
These bubbles increase the latency and reduce the throughput. The Athlon™ processor (in hard- 
ware) [28] and Itanium™ processor (in software) [5] attempt to utilize these bubbles (and increase 
throughput) by allowing other multiplications to be executed during such bubbles. 

We suggest a Booth radix-8 multiplier that allows for both operands to be either in nonredun- 
dant representation or carry-save representation. Booth multipliers with one operand in redundant 
carry-save representation were studied in [40, 8, 35]. Allowing both operands to be in redundant 
representation conceptually reduces the 3-stage pipeline to a 2-stage pipeline for all but the last 
iteration of the algorithm. 

The Booth-8 multiplier design that supports operands in redundant representation is not sym- 
metric in the sense that the first operand and second operand of the multiplier are processed dif- 
ferently during the first pipeline stage. During the first pipeline stage, operands represented as 
carry-save numbers are processed as follows: 

1. The first operand is compressed and its 3x multiple is computed. This requires two adders. 
To save hardware, we suggest employing the adder from the third pipeline stage for the 
purpose of compressing the first operand. This explains why the compression of the first 
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operand appears in Table D.2 as an operation that takes place in the third pipeline stage. 

2. The second operand can be partially compressed from carry-save representation before being 
fed to the Booth recoder. In Section D.7, we describe a recoding method that follows [35]. 

D.5.2 Directed Rounding of Carry-Save Numbers 

We propose an implementation of Algorithm 5 in which intermediate results are represented by 
carry-save encoded digit strings. Since Algorithm 5 uses directed roundings for all intermediate 
computations, we need to explain how directed rounding is applied to carry-save numbers. 

The following lemma deals with rounding up of a carry-save number. It shows that 3:2-addition 
followed by truncation after position q can be used to implement approximate rounding up. The 
3:2-addition adds a constant, called the injection, that depends only on the rounding position and 
the length of the number. 

Lemma 43 Let a[0 : t] denote a carry-save encoded digit string. Let INJCS ru(q^) G [2~ q + 
2~(9+ 1 > ) 2 • (2~ q — 2"*)] and assume that m}CS RU {q, t) is represented by a binary string I[q : t]. 
Then, 

\a[0:t}\ < \[FA(a[0:t]J[q:t])\ q \ < \a[0 : q}\ + 2~* +1 . 
Proof: We first show that 

\a[0:t)\ < \[FA(a[0:t]J[q:t])\ q \. (D.9) 

Note that from the range of INJCS bu{q, *)> li follows that I[q] = 1 and I[q + 1] = 1. Hence, 

| [FA(o[0 : t], I[q : t))\ q \ = |a[0 :q]\ + I[q] • 2"* + | [FA(a[q + 1 : *] f I[q + 1 : t])\ q \ 

= k[0:g]| + 2- q + c q -2~\ (D.10) 

where c q e {0, 1} is a carry bit to position q computed by the 3 : 2-addition. 
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We now consider two cases: 

1. c q = 1: in this case 

\[FA(*[0 : t], I[q : f])j 9 | = |a[0 : q}\ + 2"" +1 

>H0:«]|. 

2. c q = 0: since I[q + 1] = 1, a carry is not generated to position [q] only if a[q+ 1] = 0. Hence 
the \a[q + 1 : t)\ < 2~ 9 , and therefore the lower bound in Eq. D.9 holds. 

Since c q < 1, the upper bound follows directly from Eq. D.10, and the claim follows. □ 

We present in Lemma 43 the widest possible interval for miCS RU {q, t). It is probably most bene- 
ficial to use a value with the smallest number of ones, namely, iNJCSi«/(g J t) = 2~ q -f 2~ (g+1) . 

D.5.3 Optimized Implementation of Dewpoint Rounding and Back Multi- 
plication 

In this section we present an optimized implementation for computing the 2-sign (i.e. positive, 
negative, zero) of dewpoint • \B\ — \A'\. The z-sign of dewpoint • \B\ — \A'\ (combined with the 
rounding mode) determines whether the IEEE rounded quotient is RC\ or RCi. 

The following claim shows that since the dewpoint is close to the exact quotient, the absolute 
value of \B\ • dewpoint — is small. 

Claim 44 

-2" p+1 < \B\ • dewpoint - \A'\ < 2~ p+2 . 
Proof: By Assumption 40, 0 < p(N[) < 2~ p . Recall that 

dewpoint = [N WI \ p -i -f C dew (mode) 

= IK + WJrnd(mode) J p _! + Cdeu,. 
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It follows that 

dewpoint < N f k + INJ rnd (morfe) 4- Cde W (mode) 
< N' k + 2~ p+1 

The second line follows from the definition of the rounding injection lNJ rnd (mode) and the dew- 
point displacement constant Cdew(mode). The third line follows from 0 < p(N[). The upper 
bound now follows since B < 2. 

To prove the lower bound we need to show that — 2 ^| 1 + ^ < dewpoint. Since p(N[ ) < 2~ p 
and < 2, it follows that i— L < AT( + 2~p +1 . Since |S| > 1 it suffices to prove that N[ < 
dewpoint. 

We consider each of the three rounding modes. 

RZ: Let x and x 1 denote two successive representable significands that sandwich N[, namely, 
N[ e [x,x') where x 1 = x + 2~ p +\ In RZ, iNS rnd (RZ) = 0 and C^ W {RZ) = 2~ p +\ 
Hence: [N[ + lNJ rnd (raode)Jp_i = x and dewpoint = x 1 > N[. 

RI: Assume here that N[ € (x, a:'] (i.e., the equality can hold in x' instead of x). In RI, 
iNJ rnd (RI) = 2" p+1 - 2 _t and C^RNU) = 0. Hence [N( + iNJ rnd (mode)J p _i = and 
dewpoint = x 1 > N[. 

RNU: Assume here that N{ e [x + 2~ p , x' 4- 2~ p ) (i.e., we sandwich N[ between two mid-points). 
In RNU, lNJ rnd (RNU) = 2"* and C^RNU) = 2" p . Hence [ATJ + lNJ rTMf (mode)J p _i = 
x' and dewpoint = x' + 2~ p > N{. 

In all three rounding modes, dewpoint > N[ y and hence the claim follows. □ 

We use the following notation: a = dewpoint -\B\ — \A'\. Let a[— 2 : t] denote a two's complement 
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non-redundant representation of a, namely: 

t 

a = -a-2 ■ 2 2 + ^2 a i ' 2 ~ { ' 

The following claim suggests an optimized method for computing the z-sign of a. 
Claim 45 

a < 0 <^=> a p _ 2 = 1. 

a = 0 a p - 2 = 0 AND (OR(a p _ 2 , . . . , a t ) = 0), 

a > 0 a p _ 2 = 0 AND (OR(a p _ 2 , . . . , a t ) ^ 0), 

Proof: Let 

p-2 

a' = -a_ 2 • 2 2 + a i ' 2_i 
t= P -i 

By definition a" e [0, 2" p+2 ) and a' G (-oo, -2" p + 2 ] U {0} U [2" p + 2 , oo). Hence, if o! > 0, then 
a > a' > 2 _p + 2 , contradicting Claim 44. It follows that a! < 0. 

We claim that a' = 0 or a' = -2~ p+2 . Suppose that a' < -2~ p + 2 . Then a' < -2 ■ 2" p+2 . In 
this case a = a' + a ,/ <-2 - 2 _p + 2 + 2~ p+2 = -2~ p+2 , a contradiction. 

If a' = -2" p+2 , then a < 0. If a! = 0, then a > 0. Moreover, a = 0 iff o! = 0 and a" = 0. 
Similarly, a > 0 iff a! = 0 and a" > 0. □ 

We support negative numbers using a representation called two's complement carry-save repre- 
sentation [7]. In this representation, a number is represented by two binary vectors, each of which 
is a two's complement number (i.e., the most significant position has a negative weight). 

In our implementation, a two's complement carry-save representation of a is computed. For 
the purpose of computing 2-sign(a), we conclude that: (i) the carry-save digits of a in positions 
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\p — 2 : t] suffice, (ii) the subtraction of | A'\ does not require using 4 — A r as an addend; it suffices 
to use 2"fr- 3 ) - A'\p - 2 : p - 1] (see cycle 7 in Table D.2). 

D.6 An implementation with a full size multiplier 

In this section we present an implementation of our algorithm that uses a full size multiplier (i.e., 
for a target precision p, the multiplier dimensions are roughly pxp) and an initial approximation 
that whose precision is roughly p/2. To be concrete, we present detailed design for single precision 
(i.e. p = 24) that has a latency of 9 cycles. 

The main motivation for this setting is an integrated single precision and double precision 
floating-point divider. We assume that core of the design for double precision division is a half-size 
multiplier (i.e., somewhat larger than a^x (pd/2) multiplier, where p d — 53). Such a multiplier 
is larger than a 2p x p multiplier. This implies that, when performing single precision division, 
the multiplier is not only a full-size multiplier, but a double-width multiplier. We can employ a 
double-width multiplier to even reduce the latency from 9 cycles to 8 cycles (see Sec. D.6.4). 

We note that an extension of the error analysis from Section D.3 to the case of double precision 
yields the following results. Correct DEEE rounding is achieved with a "full-size" multiplier (57 x 
62 bits) and an initial approximation with precision e 0 < 2" 13 * 51 . The corresponding schedule 
for such multiplier design is listed in Table B.3. We omit the proofs and details due to space 
limitations. This double precision design achieves a latency of 11 clock cycles (including the 
initial approximation and IEEE rounding). 

D.6.1 Basic microarchitecture 

A block diagram with the four stages of our basic microarchitecture is depicted in Fig. D.l. The 
core of this microarchitecture is a radix-8 Booth recoded multiplier. The microarchitecture allows 
for feeding of previously computed products back to the multiplier as either operands. Latency is 
reduced by allowing the feedback to be in redundant representation. The multiplier also supports 
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a multiply-add operation (i.e. x • y + z). The addend (i.e. the number to be added to the product) 
is a binary number, so only one row in the adder tree is allocated for the addend. 

The multiplication circuitry is divided into 4 pipeline stages. Below we describe some details 
of the stages: 

In pipeline stage 0, an estimate for the reciprocal 1/B is computed and the values of A and B 
are compared to determine the pre-normalized value A 1 . 

In stage 1, the two operands of the multiplier (i.e., the multiplicand and multiplier) are prepared 
for the addition of the partial products in the second stage. The second operand (i.e., the multiplier) 
is recoded in Booth radix-8 digits and the partial products are generated. Following [35], the 
recoding can accept either a binary string or a carry-save encoded digit string. 

The first operand (i.e., the multiplicand) is processed as follows. The 3x multiple of the mul- 
tiplicand is computed using an adder. A feedback product, encoded as a carry-save digit string, 
can be used as a multiplicand as follows. The computation of the 3x multiple is preceded by a 
4:2-adder that computes a carry-save encoding of the 3x multiple. This carry-save encoded digit 
string is compressed to a binary number by the adder. Meanwhile, the binary representation of the 
multiplicand is computed by the adder in the third pipeline stage. This saves an adder in the first 
pipeline stage (we need to insure that the adder in the third pipeline stage is indeed available). 

In stage 2, the partial products are compressed by an adder tree. In addition to the partial 
products, there is a row that is dedicated for either (i) a round-up rounding injection (ii) an EEEE 
rounding injection (see Equation D.6), or (iii) a two's complement number (see Section D.5.3). 

Stage 3 contains an adder and a comparator. The adder has the following tasks: (i) compressing 
the multiplicands from carry-save representation to binary representation and (ii) computing RC\ 
and i?C 2 . 

The comparator is used for computing the sign of the back-multiplication. We cannot use the 
adder as a comparator due to conflicts between two pipelined divisions. 

Additional circuitry is installed in the third stage, including: (i) circuitry for feeding the adder 
with 2~ 23 and (ii) circuitry for selecting between RC\ and RC 2 . 
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D.6.2 Pipeline schedule 

Table D.2 lists the schedule of operations that implements our algorithm for single precision based 
on a "full-size" multiplier and a "half-size" initial approximation. In the error analysis presented 
in Section D.6.3 we show that a 30 x 27 rectangular multiplier together with an error bound of 
|e 0 | < 2~ 12 71 guarantees correctness. We use these parameters in the detailed description below. 
We now describe the cycles one by one. 

cycle 0: The following computations are performed in cycle 0: 

1. An initial approximation of 1/B is computed. We denote this approximate reciprocal by 
(1 — e 0 )/B. Since this reciprocal is recoded in the next cycle, it is possible to have the 
reciprocal represented in binary or redundant format. 

2. Comparison of A and B. The comparison determines A' as follows: 



A'{-\ : 23] <- i 



Ieft-shift(A[0 : 23]) if A < B 
A[0 : 23] otherwise. 



Note that since \A'\ can be in the binade [2, 4), a bit in position [—1] is required. 

cycle 1: Iteration 0 of the algorithm begins in this cycle. The approximate reciprocal (1 — e 0 ) /B 
is assigned to FLi and is fed as the second operand of the multiplier. The approximate reciprocal 
is recoded as Booth radix-8 digit string. The first operand of the multiplier D'_ x is simply assigned 
the value B. The 3x multiple of D'_ x is computed in this cycle. 

cycle 2: Two activities take place in this cycle. In the first pipeline stage, A'[—\ : 28] is assigned 
to N'_\ and the 3x multiple of NL X is computed. In the second pipeline stage the rounded up 
product DL\ - FL\ is stored in D' Q [0 : 28] as a carry-save number. The rounding up of a result 
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represented in carry-save representation is done by using the rounding up injection particular to 
carry-save representations INJCS ^(28, 28 + 14) (see Sec. D.5.2). 

cycle 3: The carry-save representation of D' 0 is compressed to binary representation in the third 
pipeline stage. In the second pipeline stage, the rounded down product NL X • F'_ x is stored in 
A^o[— 1 : 28]. Rounding down is simply achieved by truncation of the carry-save representation. 
Note that our error analysis (Claim 38) only implies that N' x < 2, hence we need position [—1] in 

Ni- 

cycle 4: In the first pipeline stage, a binary representation of F$ = 2 — 2~ 26 — D' Q is computed 
and then Booth recoded. (Recall that D f 0 is compressed to binary representation in the previous 
cycle.) In addition, the 3x multiple of Nq is computed from the carry-save representation of Nq. 
In the third pipeline stage, the carry-save representation of Nq is compressed to binary. 

cycle 5: In the second pipeline stage we compute the carry-save representation of iV INJ . 

iV INJ <— Nq - Fq + WJ rnd (mode). 

cycle 6: In the third pipeline stage, RC\ is computed by compressing the carry-save represen- 
tation of N Wi and truncating at position 23. The dewpoint equals RCi + Cdew(mode). However, 
to reduce delay, we compute the dewpoint as follows. Let carry 2 3(N MJ [24 : 28]) denote the bit 
that indicates if iV INJ [24 : 28] > 2~ 23 . A carry-save representation of the dewpoint is obtained by 
adding N INJ [0 : 23] + Cde^mode) + carry 2 3(N WJ [24 : 28]). The dewpoint is Booth recoded and 
the 3x multiple of B is computed. 

cycle 7: Back-multiplication takes place in the second pipeline stage. We use a representation 
called two's complement carry-save representation [7]. In this representation the most significant 
position has a negative weight. Our goal in the back multiplication is to compute the z-sign of 
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a = B- dewpoint — A'. (The z-sign indicates equality, greater than, or less than.) In Section D.5.3 
we show that it suffices to examine digits in positions [22 : 47] to determine the z-sign of a. 

cycle 8: In the third pipeline stage, (i) RC 2 is computed, (ii) the z-sign of the back-multiplication 
is determined from a [22 : 47], and (iii) the IEEE rounded quotient is selected from RC\ and RC 2 
according to the rounding mode and the z-sign of a. 

D.6.3 Correctness 

In this section we prove the correctness of a design based on a rectangular 30 x 27-multiplier and 
an initial approximation with a relative error |e 0 | < 2~ 12 71 . To support single precision (p = 24), 
dewpoint rounding requires that Assumption 40 holds, namely, 0 < p(N[ ) < 2" 24 . 

The analysis presented in this section is a concrete analysis for single precision. More gener- 
ally, this analysis applies to the following situation: (i) the accuracy of the initial approximation is 
roughly 2" p/2 ; (ii) the quotient is approximated by N[ (namely, the algorithm has 1 iteration of 
quadratic convergence); and (iii) the multiplier is a rectangular multiplier whose dimensions are 
not smaller than roughly p x p. 

In Table D.2 we assume that the initial approximation FL X = is represented by 15 bits. 
The following claim shows that a correct IEEE rounded quotient is computed. 

Claim 46 If the relative error in initial approximation satisfies |e 0 | < 2~ 12 n , then correct single 
precision IEEE rounding is obtained with a rectangular 30 x 27 -bit multiplier. 

Proof: All we need to show is that 0 < p{N[) < 2" 24 . To use Claim 38 we prove upper bounds 
on the relative error terms (i.e., no, ni, /o, do)- We bound these relative error terms using absolute 
error terms defined as follows: 
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Definition 14 The absolute errors of intermediate computations are defined as follows: 

nepsi = Hi ■ • F[_ x 
dep Si - di -DU FU 
fepsi ± fi*{2-D$. 

We now bound each of the relative error terms: 

Bound on do" According to Table D.2, D' 0 is represented by a carry-save encoded digit string 
with digits in positions [0 : 28], This digit string holds the rounded up value of D'^ • F!_ v This 
implies that the absolute error depso associated with the computation of D' Q is in the range [0, 2 - 
2~ 28 ) (see Lemma 43). Since the precise product D'_ x • F!_ x equals 1 — eo, this product is in the 
range 1 ± 2~ 1271 . It follows that d 0 < 2~ 27 /(l - 2~ 12 - 71 ) < 2~ 26 " 978 . 

Bound on n 0 : As in the case of D' Qj Nq is also represented by a carry-save encoded digit string 
with digits in positions [0 : 28]. This digit string holds the rounded down value of N f _ x • FLi* This 
implies that the absolute error neps Q associated with the computation of D f Q is in the range [0, 2 • 
2" 28 ). The product A! • F!_ x equals § • (1 - |e 0 |). Since A > B, it follows that A' • FL X > 1 - |e 0 |, 
and therefore, n 0 < 2" 27 /(l — |e 0 |). 

Bound on / 0 : Note that D' Q is compressed to binary non -redundant representation in the third 
cycle. We first show that fepso < 2~ 26 . The reason is that instead of computing 2 — D' 0 [0 : 28], we 
compute 2 - 2" 26 - D' 0 [Q : 26]. The absolute error is 2~ 26 - D' 0 [27 : 28]. Hence, feps 0 < 2 -26 , 
as required. Next, we expand 2 - D' 0 = 2 — (1 + d Q ) • (1 - e 0 ). We conclude that 



2 -26 

f ° - 2 - (1 + d 0 ) • (1 + |eo|) 

< 2~25.999784 
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Bound on As in the case of no, nepsi < 2 . 

We bound from below N 0 and Fq as follows: JV 0 = (1 - n 0 ) • (1 - e 0 ) • Hence, JVJ > (1 — 
n 0 ) • (1 - |eo|). Also, Fg = (1 — / 0 ) • (2 - Z? 0 ) > (1 - / 0 ) • (1 - 6 0 ). We conclude that 

2-27 

ni< (l-n 0 ) (l-|e 0 |)-(l-/o)-(l-5o) 

< 2 -26.999569 

We now apply Claim 38 to obtain: 

P(N[) < 1 - (1 - ni) • (1 - no) • (1 - fo) ■ (1 - eg - d 0 • (1 + e 0 ) 2 ) 

< 2~24.0016 < 2~24 

□ 

D.6.4 Design alternatives 

The proposed implementation is not the only way to implement our algorithm. Certain choices 
were made in order to keep the implementation simple without increasing cost by much. Below 
we list a few alternatives. 

Latency reduction. One can execute cycles 2 and 3 simultaneously and reduce the latency from 9 
cycles to 8 cycles. One way to execute the two multiplications simultaneously is to use a multiplier 
in which (i) the length of the first operand is wider than 29 + 30 bits and (ii) partial products 
corresponding to the two different products are not mixed. Note that such a multiplier is possible if 
double precision is supported (even if the multiplier is "half-size" with respect to double precision). 

The two multiplications that are executed in cycles 2 and 3 have the same second operand (i.e., 
Fij). Hence a 60 x 15 multiplier with a possibility to separate the two type of partial products 
suffices. Note that compression of D' 0 during cycle 3 is not crucial and could be skipped at the 
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expense of increasing the error term / 0 . 

Avoid the comparison A < B and using A'. In this case the exact quotient (and hence N[) 
may belong to the binade (1/2, 1). This means that N{ has to be normalized so that it is in the 
binade [1, 2). This modification has the following implications: (i) The rounding injection cannot 
be added in cycle 5 since normalization changes the value of the injection. Hence in cycle 6 an 
injection (or injection correction) needs to be added when computing RCi. (ii) The error analysis 
should be changed to deal with this case (slightly worse bounds will be obtained). 

Avoid position [—1] in the first operand of the multiplier. Position [—1] is required in the first 
operand in two places: (a) in cycle 3 to represent NL X [— 1] and (b) in cycle 5 to represent Nq[— 1]. 
If the initial reciprocal approximation has a one-sided error, namely e 0 > 0, then Nq < 2 and 
Nq[—1] = 0. As for NL X [— 1], the additional addend of the multiplier is not used in cycle 3. Hence, 
we could compute instead: mulRziN^O : 28], F'_A0 : 14]) + NL X [-1] • F1J0 : 14]. 

Save a comparator in pipeline stage 0. One could compare A < B during cycle 1 using the 
comparator in stage 3. This means that the comparator in stage 0 is not needed. Note that a 
conflict is not created between successive division operations since cycle 1 of the later operation 
corresponds to cycle 7 of the previous operation. The comparator in stage 3 is used only during 
cycle 8 so it is free to perform the comparison A < B for the later operation. 

Extra time for reciprocal approximation. In cycle 1, very little processing of FL X is required. 
This means that, if one needs more time to compute the reciprocal approximation, one can extend 
this computation so that it also takes place during part of cycle 1. 

Instant operand normalization. We chose to use a multiplier in which both the first operand 
and the product have a digit in position [-1]. This can be avoided by instantly normalizing the 
first operand before feeding it into the multiplier so that the leftmost digit is always non-zero (this 
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requires a row of MUXes). 



D.7 Implementation of the dewpoint computation 

In this section we provide additional details on the implementation of the dewpoint computation 
in the dewpoint circuit. Figure D.2 depicts the complete details of the dewpoint computation and 
recoding. 

The dewpoint is computed as the sum of 

dewpoint[0 : 24] <- iV 1NJ [0 : 23] + C<u> w {mode) + carry 2 s(N Wi [2A : 28]). 

The representation of iV INJ [0 : 28] is given as a carry-save string, and Cdew(mode) can be easily be 
generated from the input of the rounding mode as 



C few (mode) = < 



2~( p - 1 > if mode =RZ 
2~ p if mode =RNU 
0 if mode =RI. 



Note that the representation of Cdewimode) has at most one non-zero bit position that resides either 
in bit position [24] or in bit position [23]. 

The value of \dewpoint[0 : 24] | does not need to be fully compressed, but it is to be fed 
into the Booth radix-8 multiplier as the recoded operand. For the implementation we apply a 
partial compression technique from [35] that allows compression from carry-save to Booth radix-8 
recoded format through: (i.) the addition of a transformation constant (Constant = 2~ 3 + 2° + 

h 2~ 21 ) and (ii.) the application of a row of 3-bit adders with inverted most significant sum bits 

(MSSBs). 

Note that the same recoding circuitry can be used to recode the operand represented in non- 
redundant representation. We separated the circuitry for recoding only to make the explanation 



simpler. 



167 



a. 



8. 

c 



O. 



a. 
6 
o 
u 



s? 



1-1 

6 + 
i 



"S 

1 = 

a> co 
^ cs 

CM CS> 

QQ I 

^ i 

3 <m 



+ 

o 
i 



S3 -° 



+ 

CM 



3 



is 



O 



I CO 
CM 

A si 

CO 



S 

s 8 



^8 



I 



9 • 

q p 
£ a. 



168 



3 
3 



IN 
lO 
I 



ft; 



A as 
.§> g 

co <U 



ft? 8 



5 



CN 
CX 

o 



o 



o 

CO 



E 



o 

CO 



5? 



o 

co z 
o + 



So 



1 
I 

o ^ 



CM 
+ 




CO 



o 



9 

lO (J 



co bT£* 



^ o 



CO 

1 NO, 

co ^ 

*n °3 



i ssa. 



CM 

§• 8 



Is 
o o. 



169 



stage 0 



ED 



A < B? 



0 preshift 



Approx 
(1/B) 



se 



ect 



stage 1 



*5 



4:2 adder 



TT 



select 



select 3x 



3X 



rmode 



recode 
Booth-8 



dewpoint 
circuit 



select lx 



select 



stage 2 



R5 



HIE 



R6 



] E 



R7 



inject 



(10+1):2 Partial Product 
x Generation & Addition 



#cycle 



.additive input 



R8 



stage 3 



#cycle 



>-23 



0 select 1 1 select 0 



2 : 1-Add 



compute 
z-sign(a) 

rnd dec 



■ l rnd 

select I 1 



IEEE 



D' 



rounded quotient 



Figure D.l: Microarchitecture for single precision Division Implementation. 



170 







LKI i g 
[93] § f 

[«1 ^ 8 





















oo S a 
oh 

5> - ~ 

i5 S 



1^ 



8 — 



CM 

CO 



CO 
CM 



V 



Urn 

-a 

i 

CM 
CO 



5 \ 
5 \ 
5 \ 



s3S 



an8 



< 



-a 

< 



•a 
£ < 



— m 



00 .3* 

11 

CQ OS 



CO w 

p 

PQ OS 



oo _s> 

2 ° 
04 



00 

O V 

CQ Dh 



00 *2 

O oj 
CQ OS 



O QJ 

ffl as 



00 .22 

P 

03 04 



00 *i 

"is 

O u 

PQ as 



00 ^ 

11 

O u 
PQ 04 



00 .22 

P 

PQ as 



oo a> 

11 

o a> 
CQ 04 



oo « 
S w 

O u 

PQ as 



X 



I 



I 

I 



oo 



171 



D,8 Summary 

We present a complete description of an implementation of Goldschmidt's algorithm for single 
precision using a microarchitecture that contains a 30 x 27 pipelined Booth radix-8 multiplier and 
an initial approximation with a relative error bounded by 2" 12 n . This description includes detailed 
proofs of the error analysis, the rounding procedure, and various optimization techniques. 

It is possible to extend the error analysis to double precision in a setting with a 57 x 62- 
rectangular multiplier and an initial approximation with precision e 0 < 2~ 13 51 . In Table B.3, the 
schedule of such an implementation is given. This double precision implementation achieves a 
latency of 11 clock cycles (including the initial approximation and IEEE rounding). 



